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FOREWORD 
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SUMMARY 


The material presented in this report, although primarily 
concerned with the author's investigations, is representative of 
the current state of the development of nonlinear analysis tech- 
niques within the framework of the finite-element method,, Although 
the emphasis here is concerned with those nonlinearities associated 
with material behavior, a general treatment of geometric nonlinear- 
ity, alone or in combination with plasticity is included, and ap- 
plications presented for a class of problems categorized as axi- 
symmetric shells of revolution. Effects due to creep and other 
time -dependent material properties are not considered. 

The report represents an extension of two previous studies 
reported in NASA Contractor's Reports CR-803 and CR-1649. The 
emphasis of the two previous investigations was on the develop- 
ment of methods, along with sufficient computer programming and 
computation to establish their validity and effectiveness. The 
purpose of the currently reported effort is to develop a compre- 
hensive program that implements the methods and procedures of 
the two previous investigations, and thus provides practical tools 
to the designer and analyst. To accomplish this requires that the 
capabilities and capacity of the programs be sufficiently broad in 
scope so as to permit the analysis of realistic structures and 
that close attention be paid to efficiency, so that computations 
can be made more economically. 

The scope of the nonlinear analysis capabilities resulting 
from this study includes: 1) a membrane stress analysis, 

2) bending and membrane stress analysis, 3) analysis of thick 
and thin axisymmetric bodies of revolution, 4) a general three 
dimensional analysis, and 5) analysis of laminated composites „ 

Applications of the methods are made to a number of sample 
structures. Correlation with available analytic or experimental 
data range from good to excellent. 

A catalog of the finite-elements used in developing this com- 
prehensive program and a general discussion of three plasticity 
theories are presented in respective appendices. 
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1 . INTRODUCTION 


The application of finite-element methods to treat the non- 
linear behavior of structures has reached a sufficient state of 
maturity so that the results obtained from these analyses can be 
accepted with a high level of confidence. While these methods are 
generally analysis oriented, their application to industrial de- 
sign is gaining momentum in this country [1-3]. 

Although it is not customary in design practice to allow a 
structure to enter the plastic range, the need for developing a 
capability for treating the nonlinear behavior of structures never- 
theless exists. An obvious reason for this capability is that the 
prediction of failure loads (elastic and plastic instability) under 
realistic loading conditions provides a consistent, meaningful ap- 
plication of factors of safety. In addition, predicting the redis- 
tribution (whether by design or from unavoidable circumstances) of 
stresses resulting from the nonlinear behavior of a structure is a 
necessary ingredient toward its design as an efficient, safe, well 
proportioned structure. 

Compendium of Previous Studies 

While solutions exist in the literature that are applicable to 
numerous structures and materials, many of these solutions involve 
such gross simplifications that they must be used with caution. 

These simplifications were made because the complexities associated 
with plastic behavior render the solution of all but the simplest 
problems a formidable task. First, the laws of material behavior 
under complex multiaxial stress states have not been precisely de- 
fined. The stress-strain relations do not simply involve the cur- 
rent values of the stress components, as is the case for linear 
elastic behavior; rather, they depend as well upon the past his- 
tories of these components, i.e., the material has a memory asso- 
ciated with its behavior. In addition, the material laws appear 
to be different for each material. One approach to this difficulty 
is the use of approximations that provide a reasonably realistic 
representation of the essential experimentally observed features 
of inelastic material behavior. 

Another complexity associated with the general plasticity prob- 
lem is its nonlinearity. Thus, the chances of obtaining a closed- 
form solution to a specific problem are fairly remote. As was the 
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case with linear elastic stress analysis, solutions to specific 
problems had to be matched to a catalog (however limited) of pre- 
determined solutions. If the geometry and system of loads did 
not fit these known solutions, the analyst was forced to consider 
a simplified geometry and loading system that approximated the 
one in question. The nonvalidity of the principle of superposi- 
tion added an extra dimension of complexity to the nonlinear analy- 
sis. Therefore, the analyst sought numerical solution techniques 
of the repetitive type. These are usually classified as iterative 
or stepwise. In the first, the iterative approach, the entire load 
may be applied at once, and a solution is effected by using an ap- 
propriate convergent iterative technique. In the stepwise approach, 
the response of the structure is assumed to be linear over small in- 
crements of load. There are merits in both approaches, as well as 
in their combined use. 

A major difficulty with the repetitive type of solution is that 
it generally requires quite an extensive computing effort. To cir- 
cumvent: this difficulty, it is necessary not only to develop higher- 
speed machines of larger capacity but also to develop methods of 
analysis that minimize the computing effort. 

Despite these complexities, substantial progress has recently 
been made in developing general methods of plastic analysis. This 
progress is largely due to advances in the field of numerical methods 
of structural analysis, specifically, the finite-element method. The 
treatment of nonlinearities, both physical and geometric, within the 
framework of existing finite-element techniques permits analysis of 
structures of arbitrary shape and consideration of a variety of 
loading and boundary conditions. 

References 4-16 are representative of recent investigations 
concerned with incorporating the effects of plastic behavior in 
finite-element analysis. These studies describe techniques to treat 
plasticity by means of various algorithms that linearize the basi- 
cally nonlinear problem. These techniques have been divided into 
two general categories: a) the initial strain and b) the tangent- 

modulus methods. 

In the Initial-strain method, the equivalence between tempera- 
ture gradients and body forces in causing a strain field is extended 
to include plastic strains. Thus plastic effects are treated by in- 
terpreting plastic strains as initial strains. This analogy, first 
introduced In Ref„ 17 reduces the nonlinear analysis to the analysis 
of an elastic body of identical shape and boundary conditions, but 
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with an additional set of applied loads, here termed "effective 
plastic loads." The initial-strain method is advantageous because 
of the ease with which it can be implemented in a finite-element 
analysis and because the stiffness matrix developed for the elastic 
behavior of the structure remains unaffected by changes in material 
properties . 

In the tangent-modulus method (outlined in Refs. 6 and 9), the 
incremental linear constitutive relation, based on the plasticity 
theory, is introduced directly into the governing equilibrium equa- 
tions. This method requires modification of the element stiffness 
influence coefficients at each incremental load step. 

Despite differences in application of the initial-strain and 
the tangent -modulus concepts, there is, in some respects, a close 
relationship between them. References 9 and 16 discuss this rela- 
tionship and the comparative merits of the methods as applied to 
membrane stress problems, while Ref. 18 discusses the relative merits 
of each of these approaches for combined nonlinear problems. 

The development of nonlinear finite-element analysis has not 
been limited to the area of plasticity. Considerable effort has 
been directed toward the treatment of geometric nonlinearity. Ef- 
fects resulting from geometric nonlinearities alone [19-21] have 
been considered extensively, and in several instances the simul- 
taneous effects of both types of nonlinearity have been treated 
[18, 22-24]. In general, these studies have been concerned solely 
with monotonic loading conditions where the load reaches a speci- 
fied maximum value, or until structural failure occurs. Only in 
rare instances have unloading and reversed loading been considered 
[11, 16, 23]. 

In many cases, when one considers the problem of combined geo- 
metric and material nonlinearity, the solution technique used for 
the treatment of plasticity alone (i.e., an incremental approach) 
becomes especially attractive. Since the relations derived from 
the flow theory of plasticity are themselves incremental, the modi- 
fications necessary to incorporate the effects of geometric non- 
linearity are minimal. 

However, the choice of using the tangent modulus or effective 
load approach to account for material and geometric nonlinearity is 
of utmost importance because of the extensive computer times re- 
quired for repetitive solutions to these problems. In addition, 
the choice of solution algorithm (e.g., Cholesky decomposition, 
Gauss-Seidel iteration, etc.) and solution procedure (Newton-Raphson, 
first order self-correcting, etc.) in conjunction with either of the 
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above two methods can have an appreciable effect on the accuracy, 
efficiency (time), and stability of the solution. Various aspects 
of these problems have been discussed in two recent papers [18, 23], 

Grutriman- MS A Investigations 

Since 1965 NASA/Langley has partially supported several efforts 
concerning the development of methods for the inelastic analysis of 
complex aeronautical structures. A chronological summary of these 
efforts is presented in Fig. 1. The goal of these studies has been 
to provide a feasible analytical means for determining the behavior 
of structures under conditions approaching failure and for loadings 
of a realistic nature. The methods developed under these contractual 
investigations, cited in Refs. 16 and 25, are based on and represent 
an extension of, the displacement method in finite-element techniques 
of structural analysis. 

In the earliest of the previous studies (Ref. 25), consideration 
was given to the development of discrete-element methods for the 
plastic analysis of complex built-up structures in states of biaxial 
membrane stress, with particular emphasis on the effect of cyclic 
loading causing stress reversals into the plastic range. To accom- 
modate this case, the methods implemented a plasticity theory that 
can take Into account the Bauschinger effect. This theory is the 
kinematic hardening theory of Prager (Refs. 26 and 27) as modified 
by Ziegler (Ref. 28). It can represent the salient features of the 
plastic behavior of structural metals, and is readily implemented 
In a discrete-element analysis. 

The governing linear matrix equation relating nodal displace- 
ments to applied loads and initial strains was modified for the 
treatment of membrane stress states to the form of a linear matrix 
equation relating stress increments directly to increments of ap- 
plied load and initial strain. The solution to this equation was 
obtained by means of two alternative procedures. In the first, 
termed a predictor procedure, estimated values of initial (plastic) 
strain increment are used in the governing equation to yield stress 
increment. In the second, termed the stepwise linearization pro- 
cedure, the flow rule of the plasticity theory chosen is used di- 
rectly to obtain a matrix equation in terms of the unknown stress 
Increment. These two methods are discussed in detail in Ref. 4. 

The computer programs for this analysis follow a two-step procedure. 
First, an elastic finite-element program must be used to develop the 
coefficient matrices to be used in the governing equation for the 
stress increments. The program accepts nodal and member topological 
information, forms the stiffness matrix for the structure, and sub- 
sequently performs the necessary matrix operations to form the 
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coefficient matrices relating the stress increments to the incre- 
ments of applied load and initial strain. Constant stress or 
linearly varying stress elements are used in computing these co- 
efficient matrices. As a second step, computer programs are 
written, independent of the type of finite-element used, to per- 
form the plasticity analysis. 

The simplicity of this technique for membrane stress analysis 
stems from the fact that the coefficient matrices, initially de- 
veloped from an elastic analysis, remain unchanged during the course 
of the plastic analysis, and thus may be used without further 
modification . 

In the follow-on effort, reported in Ref. 16, the methods pre- 
viously developed for membrane stressed structures were extended to 
include the treatment of bending alone or in combination with mem- 
brane stresses. 

An added complexity in a plastic bending analysis is that 
plastically deforming elements cannot be treated as wholly plas- 
tic. The regions of plasticity extend through only a portion of 
the thickness, consequently, it is necessary to locate elastic- 
plastic boundaries within the thickness of the element. A technique 
developed to accomplish this both for bending alone and for bending 
in combination with membrane stresses, involves an assumed distribu- 
tion of plastic strain across the thickness (i.e., each component of 
plastic strain was assumed to vary linearly from the extreme fibers 
to an elastic-plastic boundary). A subsequent technique, developed 
under the present effort, does not require an a priori assumption 
concerning the representation of plastic strains or the elastic- 
plastic boundary through the thickness. Instead, the integrations 
necessary to obtain the components of the initial strain stiffness 
matrices are performed numerically by using the Gaussian quadrature 
scheme for in-plane integrations and Simpson's rule for the plastic 
strain integrations through the thickness. In addition, states of 
stress are evaluated at selected points through the thickness. Con- 
stitutive relations from the plasticity theory are used to evaluate 
states of plastic strain at these points, and the resulting varia- 
tion is used in evaluating the pertinent integrals associated with 
those matrices required in the plastic analysis. While the previous 
approach is quite satisfactory for cases involving monotonic loading 
conditions, it is nevertheless restrictive and leads to great com- 
plexities when consideration is given to unloading and cyclic loading 
involving reversed plastic deformation. The new technique avoids 
such complexities. 
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At the completion of Contract NAS 1-7315, the methodology had 
been extended so that the "library” of elements available for 
plastic analysis included a beam of rectangular cross section, tri- 
angular and rectangular flat plate elements, and constant strain 
and linear strain triangular membrane elements. Also included in 
this study was an investigation to determine the feasibility of in- 
cluding the treatment of geometric nonlinearity within the frame- 
work of the plasticity methods. An incremental method was used to 
account for the effects of the changing geometry of the structure 
as it deforms, and application of the combined nonlinear procedure 
was made to beam and arch structures, as reported in Refs. 13 and 16. 

The emphasis of the two NASA contracts was on the development 
of methods, along with sufficient computer programming and compu- 
tation to establish the validity and effectiveness of these methods. 
The results obtained for almost all of the sample structures con- 
sidered during the course of the contracts were compared with ex- 
perimental and/or previous analytical work, and the correlation 
ranged from good to excellent. On this basis, we believe that the 
methods and procedures developed under these contracts can be used 
with a high level of confidence. The acceptance of these plasticity 
methods as a practical means for determining the nonlinear response 
of structures subjected to realistic loadings now required that they 
be made available on a convenient basis to the structural designer 
and analyst. This requirement led to the currently reported effort. 
Contract: NAS 1-10087, the purpose of which is to develop a compre- 
hensive program that implements the methods and procedures of the 
two previous investigations, and thus provides practical tools for 
the designer and analyst. To accomplish this requires that the capa- 
bilities and capacity of the programs be sufficiently broad in scope 
so as to permit the analysis of realistic structures and that close 
attention be paid to efficiency, so that computations can be made 
more economical. 

The scope of the nonlinear analysis capabilities resulting from 
the study reported on here is presented in Fig. 2a. The capabilities 
are cataloged in terms of the types of analysis that can be treated. 
Finite-element programs, written for each of the analyses, take full 
advantage of the distinctions associated with each type of analysis 
for which they were written. The final comprehensive program con- 
sists of these programs under an executive that will selectively load 
and supervise the execution of individual programs as selected by the 
user. An organization chart of this structure is shown in Fig. 2b. A 
more detailed description of the organization of PLANS (PLastic 
ANalysis of Structures) is discussed subsequently in this report and 
in detail in a companion report. 


6 



One of the principal requirements of a comprehensive finite- 
element plastic analysis computer program is that it contains a 
library of finite-elements of sufficient variety that it is capable 
of accurately describing the state of stress and deformation exist- 
ing in a wide variety of structures. These elements should be of a 
basic shape such that arbitrary geometric configurations may be 
represented to any desired accuracy. A catalog of the specific 
finite-elements to be available in the developed plasticity program 
is presented in Appendix A. A brief discussion of each element is 
presented in this appendix in which each element is classified 
according to the type of analysis to vh ich it is applicable. 


2. GOVERNING MATRIX EQUATIONS 


The method employed in the present report uses an incremental 
formulation for the large deflection, elasto-plastic problem and is 
based on a variational principle presented in Ref. 29. The approach 
used here is identical in concept to that outlined in Ref. 24, with 
the exception that plasticity is treated by means of the initial 
strain concept [30, 31] in the present work, whereas the tangent 
modulus method [10] is used in Ref. 24. 

The use of either approach requires the use of constitutive 
relations for an elastic-plastic engineering material. In general, 
the matrix equations governing the response of a structure to some 
arbitrary history of loading are used to solve for displacement and 
total strain increment. It is necessary, therefore, to develop In- 
cremental relations linking plastic strain or stress to total strains. 
These relations are presented here in matrix form for strain hardening 
and ideally plastic behavior. A more detailed discussion is presented 
in Appendix B and Refs. 16 and 25. 


Plasticity Relations 

Appropriate incremental plasticity relations to determine values 
of stress and plastic strain developed during the histor}^ of loading 
are now considered. Hill's yield criterion [32] for an orthotropic 
material, which reduces to the Von Mises yield condition for iso- 
tropic materials, is used to predict initial yield and to obtain the 
flow rules of plasticity. The capability of handling both strain 
hardening and ideally plastic behavior is included in the analysis. 
While orthotropic behavior is included in the case of ideal plasticity, 
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only isotropic behavior is now allowed when the material strain 
hardens. There are several theories to treat the plastic behavior 
of strain, hardening orthotropic materials, but the acceptance of a 
suitable one awaits further experimental verification. A discussion 
of three theories for the treatment of strain hardening materials is 
presented in Appendix B. Two of these theories (isotropic hardening 
and kinematic hardening) have been developed for initially isotropic 
materials; the third theory (work-hardening modulii) can treat ini- 
tially anisotropic materials. 

Regardless of the theory used, linear incremental constitutive 
relations can be developed and incorporated into an equation re- 
lating generalized loads to generalized displacements of an arbi- 
trary structure. These constitutive relations can be conveniently 
represented in matrix form. 


Matrix Relations - Strain Hardening 


We can, for small strain increments, decompose the total strain 
increment {Ae T } into elastic {Ae e } and plastic {Ae} components, 

as 



Ae 


el 


+ 


1 

1 



( 1 ) 


The elastic strain increments are related to the stress increments 

{A a) by 


|Ae e | = [E ]" 1 |ao ■ 


( 2 ) 


where [E] is an array whose elements are combinations of elastic 
material constants. 


Regardless of the plasticity theory used, a linear incremental 
constitutive relation between plastic strains and stresses can be 
written. This relationship can be represented as 



[C] 



( 3 ) 


Therefore, substituting Eqs . (2) and (3) into (1) we can write 
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(4) 


'A a 


= [R ]” 1 ^Ae 1 


where [R] = [E ] ^ + [C]. 


Matrix Relations - Perfect Plasticity 


The treatment of multiaxial elastic ideally plastic behavior 
requires that the following conditions be satisfied: 

G The stress increment vector must be tangent 
to the loading surface. 


The plastic strain increment vector must be 
normal to the loading surface, where the loading 
surface is the representation in stress space of 
the initial yield function or the subsequent yield 
function after some plastic deformation has occurred. 


If f (aij) represents the yield surface, 
can be expressed analytically as 

df j „ 

^ da. . = 0 , 

da. . ij 
iJ 


the first condition 


(5) 


and provides a linear relationship among the components of stress 
increment. Thus, one of the components may be expressed in terms 
of the others. In matrix form, this can be written as 


J 

1 



[E] 



(6a) 


where [Acr] represents the independent stress components. 

The normality condition provides a linear relation among the 
various components of the plastic strain increment. This condition 
is derived from the flow rule of Eq. (B-3) , and provides a linear 
relationship in which each of the components of plastic-strain 
increment can be written in terms of any one component. This re- 
lationship may be represented in the following form 

, (6b) 

strain increment. 


9 


Ae) = [E] |a€ 


V 


where (Ae) is the independent plastic 



It is apparent from Eq. (5) that the independent increments 
of stress and plastic strain can be combined and written as the 
components of a vector, (Aco) (see Ref. 33), so that Eqs . (6) 
can be written, respectively, as 



Ae 


[E] |acd. 
[E] |aoj| 


(7a) 

(7b) 


Combining the above equations with Eqs. (1) and (2), we can form 
the following relation for the independent quantities 


Jacjd 


f T) 


= [E*] |Ae J 


where [E*] = [E ] '* 1 [E ] + [E] . 

— 


( 8 ) 


Thus, it is apparent that, upon the selection of an appro- 
priate plasticity theory to represent the nonlinear material be- 
havior, the constitutive relations can be represented in terms of 
a linear matrix relation. There are several alternatives towards 
incorporating this relation into the final equation relating loads 
to deformations. The approaches we have found to be convenient 
(in terms of required machine time) and accurate when compared with 
previous solutions and/or experiments, are outlined in the following 
subsections . 


Development of Load-Deflection Relations 

As the initial step towards the development of the governing 
matrix equation, we choose a reference state, 1^, in the body, 
for which the states of stress, strain, and deformation are known. 
We now choose the next state to be incrementally adjacent to the 
initial state with all quantities referred to the reference state, 
i.e., xi = + Au^, where xi are the new coordinates of an 

arbitrary point, are the original coordinates in the local 

coordinate system, and Auj_ are the incremental deflections of 
the point in going from the reference state to the current state 
[29]. 
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At the start of a load increment, let the stresses, surface 
tractions, and body forces acting on the structure be denoted by 
Zj_j , T®, and F?. These quantities are referred to a unit of 
"undeformed" area, i.e., before the addition of the current load 
increment. They take into account the effects of any previous ini- 
tial strains present in the body. The application of an incre- 
mental load to the body, expressed in terms of ATi and AFj_, re- 
sult in additional stresses Acrij , displacements Au-^, plastic 
(initial) strains Ae-j_j > and the distortion of the body to its 
new configuration given by x^. 

The total stresses, surface tractions, and body forces, re- 
ferred to the unit undeformed area and in the new coordinate direc- 
tions xj_, are 


a . . 
ij 

II 

h- M 

i_i. 

+ A o. 

F. 

-F° 

+ AF. 

i 

1 

i 

T. 

l 

= T° 

l 

+ AT. 

i 


The development of the governing matrix equation may be ap- 
proached by one of several alternative procedures . The authors 
choose here the principle of virtual work, which, for an incre- 
mental method, may be written as [29]: 


p 


( 2. . + Aa . . ) 6 (Ae . . ) dV 
lj ij 

^ V 


(T, + AT i )6(Au i )dS 




(F° + AF i )5(Au i )dV 


'V 


( 10 ) 


Here Ae^j is Green's strain tensor that refers to the original 
or "undeformed" volume of the element 


Ae . . = Ae . . + An . . 
1 J 


( 11 ) 
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In this expression Ae-y are the terms that represent a linear 
strain-displacement relationship, while A-q -j_ j are those associated 
with the nonlinear terms in the strain-displacement relationship* 
The incremental constitutive equations are taken to be in the fol- 
lowing form 



E ijkl 


(A£ kl 


Ae kl } 


( 12 ) 


where Ae^ are the initial or plastic strains developed in the 
current increment based upon the "undeformed” geometry. These 
are assumed to be small and independent of the total strains. The 
terms, are the linearly elastic material properties. 

Substituting the stress-strain relations, Eqs . (12) and (11) 
into Eq. (10) yields 


l Ae ij E ijki 5 (Ae ki ) + V s ^ij)]^ 


V 


AT . 5 ( Au . ) dS + 

l v l 


AF . 5 (Au . ) dV 

l 1 


V 


Ae . .E. ., , 6 (Ae. AdV - 
xj ijkl v kl 


'V 


2. . 6 (Ae . . ) dV 


V v 


(13) 


■> 

\ 

T?6 (Au^) dS - 

F?6 (Au ± ) dV 

J s 

V 


A ^ij E ijkl 6(Ae kl )dV + 


V 


Act. . 6 (At) . . )dV 
ij ‘xj' 


V 
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We now have an equation that is similar in form to that pre- 
sented in Ref. 24, with the exception of those terms associated 
with initial strains. As in Ref. 24, it is assumed that, although 
total strains may be large, incremental strains from one neighbor- 
ing state to the next are small, and thus the last two terms of 
Eq. (13) (which are cubic and quartic in displacement increments) 
may be neglected when compared to terms that are quadratic in dis- 
placement increments. These terms that are neglected lead to the 
matrices [Nq] and [N 2 ] of Ref, 34 and must be retained in a 
total Lagrangian formulation. An additional matrix due to the 
presence of initial strains is also generated from the last term 
of Eq. (13), but as it contributes terms of the same order of mag- 
nitude as the other termsneglected, it too need not be retained. 

We then have 


■a 

[ 2 . .6 (An . .) + Ae, ,E. ., ,6(Ae. . )ldV = 
ij 'lj' kl ijkl ij / J 

J V 


AT. 6 (Au . ) d3 + 
x 1 


AF (Au_. ) dV 


V 


+ 


Ae . . E . . , , 8 ( Ae. . ) dV 
xj ijkl x kl' 


V 




■5 



T?5 (Au.. ) dS + 

F?5 (Atm) dV - 

2. . 8 (Ae . . ) dV 

V 

S 

V J 

V 


(14) 


If the initial stress state, 2qj , Tq, and Fq, is in equilibrium 
at the start of the incremental step, then the last three terms of 
Eq. (14) vanish and we get the following incremental initial strain, 
large deflection relation: 
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j 2 . .5 (An . . ) + Ae. E . . , 6 (Ae . . ) ] dV = 

[ ij ‘ij 7 kl xjkl xj \ 


V 


AT i 6(Au i )dS 


(15) 


+ 


AF jL 5(Au i )dV + 


V 


Ae . .E. ., , 6(Ae. n )dV 
ij ijkl v kl 7 


V 


The first term of this equation yields the initial stress 
stiffness matrix after the rotations (or other nonlinear terms) 
have been expressed in terms of nodal degrees of freedom. The 
second term leads to the conventional stiffness matrix. The first 
two terms on the right side lead to the consistent load vectors 
for surface tractions and body forces, respectively.. The last 
term on the right side leads to the initial strain stiffness matrix 
which is multiplied by a vector of plastic (initial) strains to be 
used as an "effective" plastic load vector. 

Because we use a predictor procedure, however, the initial 
stress state may not be in equilibrium before the current load step. 
The results for the next step may be adjusted or corrected for this 
imbalance by introducing a residual force given by [24]. 



T?6 (Au £ ) dS + 


F i B ^ Au i^ dV 


o 


V 


2. . 5 ( Ae . . ) dV 
ij 


V 


(16) 


Any discrepancies due to the neglect of the change in direc- 
tion of the load are also accounted for in Eq. (16), because the 
total load is applied to the structure in its current configura- 
tion. The total stresses o-h obtained at the end of load incre- 
ment N become the initial stresses for step (N + 1) . These 
must now be related to the new deformed area (which is the unde- 
formed area for step. N + 1) . The transformation that accomplishes 
this is presented in Ref. 24 and written here as 


2 . . 

Kl 


Ae, , ) o . . 
kk ij 


+ (ie jk + to Jk ) 0 . k +(ae. k + /to lk )a. k 


(17) 
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where the Acdj^ are the incremental rotations. Similar trans- 
formations must be carried out for the surface tractions, body 

forces, and initial strains. 

We will, at this point, mention that the last two terms of 
Eq. (13) need not be neglected and can be included without the 
formation of any additional stiffness matrices besides the re- 
quired initial stress stiffness matrix and initial strain sitff- 
ness matrix. These terms may be retained in a predictor process 
in which values for Aaij and Arj j obtained from the previous 
step (appropriately extrapolated) are used in the formulation of 
the appropriate stiffness matrices for the next incremental step 
This should permit the use of larger step sizes in the current 
formulation. This latter concept was not used to obtain the re- 
sults presented here. 


Final Matrix Equations and Solution Procedures 


The final form of the incremental equations used in the finite 
element formulation is obtained from Eq. (14). The displacements 
(Au), linear total strains {Ae-yl, and rotations { Acoj[_ j } are re 
lated to the nodal generalized displacements {Au^} via the follow 
ing matrix relations: 


{ Au } = [ N ] { Au^ } 

{Ae } = [W]{A Ui ) 
(Acp..) = [f2]{Au. } 

i-J 1 


We then have 

([k°] + [k 1 ])fAu.) 


{ AP^} 4- {AQ^} + {R_^} 


where 


[k°] 


[W] i [E] [W]dV 


( 18 ) 


(19) 


V 



[k 1 ] = [n] T [2 ] [n]dv 

J v 


(AP i ) = [N] T { AT°]dS 

’ S 


{ AQ ± } = [W ] T [E ] { Ae }dv 

^ V 


{R i } = [N ] T { T°} dS - 

and body force terms have been neglected, as they are not considered 
in this report. Here [k^] is the conventional stiffness matrix, 

[k 1 ] is the "initial stress" or geometric stiffness matrix, (AP^} 
is the vector of applied loads, fAQ.j_} is the effective plastic 
load vector, and [R-jJ is the vector of residual forces due to the 
existence of any equilibrium imbalance that may exist because of 
the predictor nature of the numerical solution procedure. 

Equation (19) is derived by using a moving coordinate system 
fixed to the body. It is valid for large elastic-plastic deforma- 
tions, provided the appropriate nonlinear terms are retained in the 
strain-displacement relations and the total strain increment can be 
simply decomposed into elastic and plastic components. Additionally, 
proper transformations from the previous to current coordinate sys- 
tems must be used so that the changes in orientation and volume of 
the elements must be accounted for [23, 24]. In this report, con- 
sideration is restricted to small strain, moderate rotation prob- 
lems. Alternative solution techniques to solve Eq. (19) are pre- 
sented and comprehensively discussed in Ref. 33. For completeness, 
a brief summary of these techniques is presented here. 
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Material Nonlinearity 

If we initially neglect the effects of changing geometry and 
the equilibrium correction term, and write the elastic stiffness 
matrix as 


[k°] + [k 1 ] = [k ] , ( 20 ) 

and use the superscript "k" ' to denote the "k th " load increment, 
then Eq. (19) may be written as 

[k g ] k {A Ui } k = + (AQ i } k . (21) 

It should be noted that the initial stress stiffness matrix [k 1 ] 
is not necessarily associated only with geometric nonlinearity; it 
may also be required in such other cases as the bending of plates 
subjected to membrane stress. 

In Ref. 33, Eq. (21) is put into alternative forms suitable 
for numerical solution; distinctions among the basic forms employed 
are accomplished by using the term ''method," and distinctions among 
the solution procedures are effected by using the term "procedure." 

Displacement Method — Predictor Procedure . In the first 
"method" to be treated, the "effective plastic loading" is taken to 
be equal to that computed in the preceding load increment, and is 
thus taken as a known quantity In the equation. 

The final form of Eq. (21) is 

[k e ] k [& u .} k = {AP.p + faQ.p' 1 , (22) 

where k-1 is the preceding load step. The use of this type of pre- 
dictor procedure obviates ’the necessity of introducing the plastic 
stress-strain relations explicitly into the governing matrix equation. 

The incremental solution technique using Eq. (22) reduces to a 
sequence of linear problems in which the applied loading Is con- 
stantly modified by the effective plastic load vector. Thus, with 
the increments of generalized displacement obtained from Eq , (22), 
the strain-displacement relations and Eqs . (4) and (3) are used to 
obtain the complete solution for increments of total strain, stress, 
and plastic strain, respectively, assuming elastic strain-hardening 
material behavior. The corresponding relations [replacing Eqs. (4) 
and (3)] for an elastic, ideally-plastic material are given in 
Eqs. (7) and (8). After summing all incremental quantities to de- 
termine current values of the pertinent variables, new values of 
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the increments of fictitious load [AQ] are determined for each ele- 
ment in the plastic range, and the procedure is repeated until the end 
of the loading process is reached, A further discussion of the use of 
the effective plastic load concept in conjunction with the basic gov- 
erning matrix equation, Eq. (22), can be found in Refs. 13, 16, and 25. 

Strain Method — Predictor Procedure . The predictor procedure 
solution technique can also be applied in an alternative formulation 
of the problem involving a direct solution for the increments of total 
strain. This alternative formulation is applicable to those problems 
in which an explicit solution for displacements, or their increments, 
is not required, and where the initial strain stiffness matrix [k*] 
does not change throughout the loading range. The governing matrix 
equation in this formulation is determined by substituting Eq. (20) 
into Eq . (19), and making use of the nodal strain-displacement rela- 
tion of the following form: 

{AeT} = [A] { AP. } + [J] { Ae . } , (23) 

r i 

where 

[a] = [w Kk r 1 , 

j_ C 

and 

[J] *= [W.][kJ [k ] , 

j_ t: 

•with 

[k*] - f [W] T [E]dV 

J y 

where [W^] represents the matrix relating the linear component of 
the strain-displacement equations to the nodal degrees of freedom. 

If we wish to use a predictor procedure to solve Eq. (23), we 
must write this relation in the following form: 

X 

[AeT] = [A] { AP . } 1 + [J]fAe J 1 " 1 . (24) 

Using values of { As } estimated in this way, we can find the un- 
known total strain increments from Eq. (24), and then find the in- 
crements of stress and plastic strain from Eqs . (4) and (3) for strain- 
hardening behavior, or from Eqs . (7) and (8) for ideally plastic be- 
havior. In Refs. 4, 6, 7, and 25 the predictor procedure is also 
formulated in terms of a governing matrix equation relating increments 
of stress to increments of lead and plastic strain, similar to that of 
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Eq. (23). This procedure, referred to as a "constant stress" pro- 
cedure, has beer ohowr to lead to a characteristic numerical in- 
stability [4, 7, z5] . No such instability occurs when the pre- 
dictor procedure is used in conjunction with Eq. (22) or Eq. (24). 

The predictor procedure, involving the use of estimated values 
of plastic strain in Eq. (22) or Eq. (24), has computational ad™ 
vantages since the solution requires only matrix multiplication in 
each load step once the corresponding effective plastic load vector 
is formed, provided the matrix [k e ] is constant and thus need be 
inverted only once. This differs from r.he direct substitution pro- 
cedure, to be discussed below, in which matrix inversion or simul- 
taneous equation solution is required at each load step. However, 
a disadvantage associated with the predictor procedure solution 
technique is a "drifting" of the numerical results from the true 
solution as plastic strain proceeds. One may choose to use small, 
load increments for improved accuracy, thereby reducing the compu- 
tational advantage of this procedure, or the results for any incre- 
ment may be adjusted or corrected toward the true solution by intro- 
ducing a residual force {R-p3 to ensure equilibrium of the total, 
system. At any step, the residual force may be computed from the 
governing equations written in terms of total quantities. Thus, 
for the displacement method, we have 

C R - } k = - [k ] k fu } k + [P.} k + (Q. } k “ 1 , (25) 

JL t- JL JL X 

and for the strain method 

k • 

(R.] k = - [eh + [A] {P.} k + [J] U .} 1 "" 1 • ( 26 ) 

i. JL JL J- 

The value of the residual force is used in the next incremental 
load step, i.e . , 


[k ] k (Au } k = (AP } k + [AQ } k ” i + {R } k . (27) 

e 1 1 1 1 


Displacement and Strain Method — Direct Substitution Procedure . 
The use of Eqs , (22) or (24) is usually associated with the initial 
strain method of finite-element plasticity analysis. An alternative 
approach, commonly referred to as the tangent modulus method, in- 
volves the direct substitution of the incremental constitutive plas- 
ticity relations into the governing matrix equation, Eq. (21). For 
an elastic, strain-hardening material, Eqs. (3) and (4) may be 
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combined to yield an incremental relation between plastic strains 
and total strains. For an elastic, ideally-plastic material the 
incremental relation between plastic strains and total strains 
is obtained from Eq. (8). Thus, regardless of material behavior, 
an equation of the following form may be written at those nodes 
or elements in the plastic range: 

(Ae^ = [S] {AeJ} , (28) 

where [S] reflects strain-hardening or ideally-plastic behavior 
and is a null matrix at those nodes or elements in the elastic 

range . 


An incremental plastic strain-displacement relation is ob- 
tained by substituting the total strain-displacement relation 
into Eq. (28), so that we can write 

{Ae.} = [S][W i ] {Au i } . (29) 

A linear incremental load-displacement relation of the form 
[AP i ) = [k^J [Au^] , (30) 

where [top] = [k e ] - [k ][S][Wj_], results from substituting 
Eq. (29) into Eq. (19). Alternatively, after extensive manipula- 
tion (see [33]), 

[kj,] = [kp] + [k 1 ] , (31) 

where 

[k ] = fff [W.]'[M][W.]dv , 

^ v 

[M] = [R] ^ for strain-hardening behavior 

r * “1 

= [E J for perfectly plastic behavior. 


A direct formulation of the equations in the tangent modulus form 
has been presented in [22, 24]. 



The matrix [kp] in Eq. (31) may be regarded as a "plastic 
stiffness matrix" since it explicitly contains the effect of plas- 
ticity and enters into the analysis as an additional component of 
the total stiffness matrix. Further, since the elements of [k 1 
are functions of the instantaneous stress state, they must be 
evaluated at each incremental step. Similarly, substituting the 
incremental plastic strain-total strain relation of Eq. (29) into 
Eq. (23) results in the following equation: 

{AeJ} = [Y] {AP_ l } , (32) 

where [Y] = [A] + [J][S]. 

The direct substitution procedures, as is the case with the 
predictor procedures, may make use of an equilibrium correction to 
account for errors associated with the linear incrementation of 
the problem. For example, in the displacement method we can write 

[kp] {AUp} = { APp) + [Rp] , (33) 

where [Rp] = [kp] [up] - (Pp). 


Summary of Methods for Plastic Analysis 

A distinction among the various formulations is associated with 
the solution procedures used, which may be named the predictor and 
the direct substitution procedures. In the former, estimated values 
of plastic strain are used in the governing linear matrix equation. 
Thus, plastic effects are treated in the linear matrix equation by a 
modification that is external to the stiffness influence coefficient 
matrix. In the direct substitution procedure, plasticity is accounted 
for by means of an "internal" modification of the stiffness matrix. 

The direct substitution or internal modification procedure, 
though it retains the errors associated with stepwise linearization, 
does eliminate the propagation of error associated with the pre- 
dictor or external modification procedure. This improvement in 
accuracy for a given magnitude of the load increment is, however, 
accompanied by an increase in the number of numerical operations 
required to obtain a solution. These operations can be computa- 
tionally expensive, since the elements of the influence coefficient 
matrices, [kp] of Eqs. (30) or (31) or [Y] of Eq. (32), must be 
recomputed at each incremental step of loading. The effect of this 
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can be mitigated by increasing the magnitude of the load increment, 
but at the cost of greater inaccuracy. A choice between the two 
basic procedures thus involves a tradeoff between smaller load in- 
crement but less computation per increment, in the case of the pre- 
dictor procedure, and larger load increment but more computation per 
increment, in the case of the direct substitution procedure. This 
choice will not be obvious in any given problem. 

An approach that combines the two procedures might prove to 
be the most effective. For example, the predictor procedure may 
be sufficiently accurate in those regions of a structure where 
plastic flow has begun but has not yet been substantially developed. 
In those regions where plastic effects are predominant, the direct 
substitution procedure could be used. 


Implementation of Combined Plasticity and Geometric Nonlinearity 

A primary consideration in choosing a method for the analysis 
of geometric nonlinearity from among the several currently avail- 
able is the ease with which it can be combined with methods of 
plastic analysis. For this reason, our approach is based upon a 
linearized incremental formulation, i.e., one in which the non- 
linear analysis is reduced to the solution of a sequence of linear 
incremental equations. In Refs. 6, 10, and 35, this approach was 
used to solve problems involving geometric nonlinearity. Since the 
plasticity relations are themselves incremental, and the methods 
discussed in the preceding section depend upon a revision of the 
governing matrix equation in each loading step, the modifications 
necessary to incorporate "large deflection" terms are minimal. 

The method of solution of the small strain, geometrically non- 
linear problem discussed here involves the solution of a sequence 
of "beam-column" type problems, using Eq. (19), in which values of 
the membrane stress resultants and the geometry of the deformed 
structure are updated in each increment of loading. For sufficiently 
small loading increments, the increments of rotation in any finite- 
element will be small as measured with respect to a local coordinate 
system that translates and rotates with the element in successive 
loading steps (but is assumed to remain fixed within any one loading 
step) . Consequently, squares and products of the increments of 
rotation may be neglected in computing increments of membrane strain 
[29]. 


Because of the presence of geometric nonlinearity, the entire 
element stiffness matrix [k e ] in Eq. (20) must be reformed in 



each loading step, with current stress levels and geometry being 
used. In the discussion of the development of the elastic stiff- 
ness matrix, it was mentioned that the only component matrices re- 
quired are the conventional stiffness matrices (those not dependent 
upon the presence of stress) and the initial stress stiffness matrix 
The latter matrix accounts for the change in bending stiffness due 
to the presence of membrane loads. In the development of the ini- 
tial stress stiffness matrix, the membrane stresses taken into 
account are those present at the beginning of the loading step, any 
further changes in these stresses occurring during the loading step 
being neglected in that development. This constitutes the lineariza 
tion of the procedure during an increment of loading. 

Some investigators (see [21, 22]) have indicated the need for 
an additional matrix, termed "the initial displacement" matrix, for 
the treatment of geometric nonlinearity. Because the current analy- 
sis uses a "moving" local coordinate system, this additional stiff- 
ness matrix is not required (see [24] for a more complete discussion 
of this point) . 

For the large deflection elastic-plastic analysis, the load is 
applied in small increments from the initial unloaded state. At 
the end of each increment, new increments of deflection, stress, 
strain, and plastic strain are calculated. Total quantities, such 
as the initial stresses, are calculated by using appropriate trans- 
formations, and the geometry of the structure is updated. Again 
the plastic strain increments used are those calculated in the pre- 
vious step. The total stiffness matrix is reformed at every in- 
crement, together with the incremental load vector, plastic load 
vector, and residual load vector. The element contributions are 
then assembled and the system of linear incremental equations is 
again solved and the process repeated until the maximum specified 
load is reached or structural failure occurs. If the response to 
cyclic loads is desired, the load increment is reversed at the maxi- 
mum load, and the incremental process is repeated until the new 
maximum (minimum) load is reached. This procedure is then repeated 
for as many load cycles as desired. 

For the large deflection problem, the most time-consuming fea- 
ture is the reassembly of the stiffness matrix and solution of the 
linear incremental equations. It becomes convenient, therefore, to 
consider the possibility of treating the large deflection terms as 
well as the plasticity effects as effective loads. This may be done 
by rewriting Eq. (19) as 
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[k°] {A Ui } k = - [k 1 ] {Au i } k_1 + (AP i ) k + (AQ i } k " 1 + {R.} , (34) 

where the product of the initial stress stiffness matrix and the 
vector of displacement increments of the previous step is now 
treated as an "effective geometric load." The stiffness matrix 
[k°] may be re-formed every M steps (M > 1), with the possi- 
bility of saving considerable time. The use of this solution pro- 
cedure may lead to numerical instabilities when the nonlinearities 
become large [20] s although none were observed in the limited num- 
ber of problems solved by the authors. The use of the geometric 
terms as effective loads is not new, and has been used in many 
Lagrangian formulations with great success [20] . 

If large deflection terms are not important, considerable sim>~ 
plification results. For the problem of material nonlinearity 
alone If we consider the displacement method-predictor procedure, 
the sequence of computations followed is to calculate the value of 
the load at which plastic deformation first occurs. From this 
point, the load is incremented to a maximum value, with new in- 
crements of displacement, plastic strain, and stress calculated at 
each step and total values obtained by summing incremental values. 
The plastic strain increments used in the plastic load vector are 
those calculated from the previous step. Because this is a small 
deflection analysis, the stiffness matrix need never be re-formed. 
The residual force vector was not used In any of the small deflec- 
tion problems presented although it could have been so used. At 
the maximum load, a new critical load for which yielding begins in 
the reverse direction is calculated, based upon elastic unloading 
to this point. Procedures for determining this load are presented 
In Refs. 16 and 25. This critical load may occur before all the 
load is removed from the structure because of the presence of 
residual stress and the existence of the Bauschinger effect. At 
this new critical value, the load is incremented to the new speci- 
fied maximum (minimum) value, and this procedure is repeated for 
as many half cycles as desired. 


3. APPLICATIONS OF NONLINEAR METHODS 


The methods and procedures of the previous section have been 
applied, to a variety of problems during the course of several years 
of investigation by the authors. Comparison between the predictor 
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procedure and the direct substitution procedure involving plasticity 
alone is made and evaluated for the strain method as applied to mem- 
brane stress problems [25], and for the displacement method for 
problems involving bending stresses [16]. A survey and thorough 
evaluation of procedures for combined geometric and material non- 
linearity is presented in Ref. 36. 

The development of practical methods of treating the nonlinear 
response of structures within the framework of finite-element tech- 
niques has required that they be made available on a convenient 
basis to the structural designer and analyst. This requirement 
has led to the development of a comprehensive program that provides 
a practical tool for the designer and analyst. A brief description 
of this program is presented in the next section and, in more de- 
tail, in a companion report. To accomplish the development of such 
a program requires that the capabilities and capacity be sufficiently 
broad in scope to permit the analysis of realistic structures and 
that close attention be paid to efficiency, so that computations 
can be made economically. 

On the basis of the authors' experience we have found that for 
plasticity problems alone, the displacement method-predictor pro- 
cedure is the most attractive in terms of generality and ease of ap- 
plication. Furthermore, it is highly competitive in terms of effi- 
ciency and accuracy as compared with various other algorithms dis- 
cussed in this report. Hence, in the program developed in this study 
and for the results to be subsequently presented, this technique has 
been used exclusively. For problems involving combined material 
and geometric nonlinearity the predictor procedure is used for both 
geometric and material nonlinearity until the geometric nonlineari- 
ties become "large." From then on the geometric nonlinearities are 
included via the tangent modulus approach while the material non- 
linearities are still handled via the effective load technique. 

The equilibrium correction term is used in each increment. This 
procedure is effective [23] although not necessarily the most ef- 
ficient [18]. 


Elastic-Plastic Analysis of Representative Two and Three Dimensional 
Structures 


Our basic philosophy of program organization for material non- 
linearity alone, as presented in the next section, is to separate 
classes of analysis into individual groups, with each group defined 
by the physical problem to be solved. For example, groups can be 
defined for membrane structures, combined bending and membrane 
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structures, thick or thin bodies of revolution, general three 
dimensional solids, or laminated composites. On the basis of 
this concept, each group is in itself an independent finite “element 
computer program, with an associated element library that can be 
individually loaded and used to solve the problem class of in- 
terest. In this manner any simplification or specialization ger- 
mane to the individual analysis can be incorporated. 

Each group can contain subprograms common to all other pro- 
gram groups as well as subprograms distinct to its group alone. 
Linkage to any individual program group is made through an execu- 
tive program that performs the function of loading the individual 
analysis programs. Data set management and core allocation are 
handled by each of the group programs, and linkage to any group 
can be made only through the executive routine. 

Representative problems from each of the individual analysis 
programs are presented to demonstrate the scope, accuracy, and 
general adaptability of the analysis techniques. 


Membrane Stress Analysis 

The element library for this program contains the following 
elements: 1) constant strain triangle, 2) linear strain triangle, 

3) hybrid triangles linking the constant strain and linear strain 
triangles, 4) warped quadrilateral shear panel [37], 5) stringer 

elements, and 6) a beam element. This program Is designed to be 
used for the analysis of arbitrarily shaped structures (stiffened 
or unstiffened) in which the loads are those that generate pri- 
marily membrane stresses. 

Uniformly Loaded Sheet with a Central Crack . An illustration 
of the application of the membrane stress program is that of a uni- 
formly loaded sheet with a central crack. The idealization of a 
quadrant of the sheet, as shown in Fig. 3, involves 455 elements, 
219 vertex nodes, and 250 mid-side nodes, resulting in 902 
degrees - of -'freedom when symmetry boundary conditions are applied. 

A mixture of linear strain elements (used in the region surrounding 
the crack tip), transition elements, and constant strain elements 
are used to make up the network shown. 

Results from the analysis in the form of the distribution 
along the horizontal axis of symmetry of the stress component in 
the direction of loading is shown in Fig. 4a. Results from the 
finite-element analysis for the elastic behavior compare favorably 
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with a continuum solution [38]. Results for the plastic behavior 
are shown for two levels of loading corresponding to stress levels 
for which the gross section stress (a G ) is 37.6 percent and 
75.3 percent of the yield stress of the material. The effect: 
of plasticity in decreasing the stress gradient in the vicinity 
of the crack tip is clearly evident from the figure. Also shown 
are the residual stress distributions that result from complete 
unloading from the respective maximum loads of 37.6 percent and 
75.3 percent of yield stress. The residual stresses in the 
vicinity of the crack tip are significant and can be expected to 
represent an important factor in determining the fatigue life of 
such structures. 

The displacement profile of one-half of the crack is shown in 
Fig. 4b. Once again, a comparison of the finite-element results, 
assuming elastic behavior compares favorably with the continuum 
solution. The results from the finite-element analysis for load- 
ing in the plastic range indicate that the displacement profile is 
virtually unchanged from the elastic profile. This is attributed 
to the fact that plasticity is localized to a rather small region 
around the crack tip. That is, for monotonically increasing loads, 
the localized region of plasticity does not appreciably affect the 
deflection profile to the extent that it influences the stress dis- 
tribution in the vicinity of the crack tip. However, because of 
the development of a plastic zone there is a residual displacement 
profile along the entire crack length as shown in Fig. 4b. This 
profile has a "dog-bone" type pattern, i.e., the maximum opening 
occurs near the crack tip and becomes uniform toward the center of 
the crack. 

The propagation of the plastic zone in the uniformly loaded 
sheet with a central crack is shown in Fig. 5a for monotonically 
increasing loading and in Fig. 5b for unloading from the maximum 
load. The zones are determined by drawing a curve through the 
centroid of those elements that are plastic at the corresponding 
load. As shown in Fig. 5b, prior to the removal of all the load 
(oo/cTyieTci = 0.345) a compressive plastic zone develops in the 
region of the crack tip. The plastic zone that remains upon the 
complete removal of the loads is also shown in this figure. 

Additional problems, typical of the type that can be treated 
by this module of the program, have been presented in previous 
publications by the authors (i.e., 11, 16, and 25). Although not 

illustrated here, the current capability includes tne treatment 
of general out-of-plane stiffened and unstiffened structures. 
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Bending and Membrane Stress Analysis 

This module of the comprehensive program is to be used for 
the general analysis of structures in which bending and membrane 
stresses are generated. Its finite element library consists of 
flat triangular bending [39] and membrane elements, and a general 
beam element to represent stiffeners. 

A wide range of applications of prototypes of this program 
appear in several of the author's previous publications (i.e., 

13, 16, and 40). The problems considered to date include ini- 
tially flat rectangular, circular, and annular plates subjected 
to loads causing bending alone and combined bending and membrane 
stress states. The current capability includes the treatment of 
general out-of-plane stiffened and unstiffened structures. 


Axisymmetric Analysis of Bodies of Revolution 

The element library for this program consists of a revolved 
triangle [41] for general thick walled bodies of revolution, an 
isoparametric shell element [23, 42], and a ring element [43]. 


Full Sphere . The problem of a full sphere subjected to a 
"cos 2 e" loading is used to illustrate the accuracy of the elastic 
solution obtained by using the revolved triangles or the shell ele- 
ment. An idealization of a quadrant of the sphere, as shown in 
Fig. 6, consists of 153 triangular elements, 95 nodes, resulting 
In 183 degrees-of -freedom when symmetry boundary conditions are 
applied. An idealization using the shell element consists of 
15 elements, 16 nodes and results in 60 degrees-of -freedom. The 
meridional segments shown in Fig. 6 represent the meridional lengths 
of the shell elements. 

A comparison of results for the radial displacement distribu- 
tion versus the meridional direction is shown in Fig. 7. The re- 
sults from the two separate finite-element solutions are compared 
with an elasticity solution [44]. The correlation of results be- 
tween the three solutions is quite good, with a maximum divergence 
occurring at, and in the neighborhood of, the apex. A similar 
correlation exists when a comparison is made for the circumferen- 
tial stress versus the meridional direction, as shown in Fig. 8. 
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Thick Tube Under Uniform Internal Pressure . The elastic - 
ideally plastic behavior of an internally loaded circular tube was 
considered to be a representative test of the accuracy of the non- 
linear analysis using the triangular elements. Results from the 
finite-element analysis are compared to the results from a solu- 
tion obtained in Ref. 45, and shown in Figs. 9-11. For this par- 
ticular case, the ratio of the outer to inner radius equals 2. 

The external radial displacement and radius of the elastic-plastic 
boundary versus the internal pressure p are plotted in Fig. 9. 

The comparison of results between the two solutions is good for 
the entire range of loading considered. The resulting radial and 
circumferential stress distributions for various positions of the 
elastic-plastic boundary are shown in Figs. 10 and 11. The corre- 
lation with the results of Ref. 45 is again seen to be quite good, 
despite a discontinuous distribution of circumferential stress at 
the elastic-plastic boundary as shown in Fig. 11. It should be 
noted that the stresses plotted represent the "average" stress at 
a node; that is the average value of the stresses in all those 
elements common to a node. 

Sheet with an Oversized Fastener . Another illustration of the 
use of this module of the program considers the residual stress 
field surrounding a hole in a sheet into which an oversized fastener 
has been driven. As shown in Fig 12, the "head" of the fastener, 
which must be driven through the hole, is tapered to a maximum 
diameter Dp, with the shank stepped down to a diameter, D 2 . Re- 
sults for the residual circumferential and radial stress distribu- 
tion along a radial line on the outer surfaces of the sheet are 
shown in the figure for three cases of prescribed interference, 

Dl -2a. These stresses apply with the fastener shank in the hole. 
Linear strain hardening material behavior is assumed for the sheet 
material in the plastic range. It is significant to note the small 
increase in favorable compressive residual stress at the edge of 
the hole as compared to the substantial increase in tensile stresses 
at some distance from the edge. For the case involving the largest 
interference, the level of this tensile stress was found to be un- 
acceptable from both fatigue and stress corrosion considerations. 

Torispherical Shell . The accuracy of the shell element, used 
in this module of the PLANS program, is demonstrated in Ref. 42 in 
the form of an extensive number of applications and comparisons 
with existing solutions. An application of the use of this element 
for a plastic analysis is presented in Fig. 13 where the load versus 
apex deflection for a torispherical shell under uniform Internal 
pressure is presented for various load increments, and a comparison 
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Is made with the results obtained in Ref. 46 where elastic-perfectly 
plastic behavior was assumed. The results from Ref. 46 were obtained 
by using; load increments of 1.03 x 10 4 N/m , and are virtually 
identical with those of the present analysis, where a load incre- 
ment of 2.76 x 10 2 N/m 2 was used. As seen in the figure, halving 
this load increment produces a significant change in the results 
only at a load above the theoretical collapse load predicted by 
limit analysis [47]. The use of the initial strain method, wherein 
the plastic behavior is accounted for by an "effective plastic 
load” vector, requires smaller load increments than a tangent 
modulus method [16]. However, increment size alone is not the 
sole criterion governing the efficiency of one method versus 
another. The increase in computing time associated with the use 
of smaller increments in the initial strain method is offset by 
the fact that the stiffness matrix need never be reformed after 
the first step. Additional evidence that indicates that the ini- 
tial strain procedure is competitive from the standpoint of com- 
puter time requirements is presented in Ref. 48. 


Figure 14 shows the load versus apex deflection curve for the 
same shell for one full cycle of loading. The load is varied be- 
tween amplitudes of ±80 psi (±5.52 x 10^ N/m 2 ) and back to zero. 
It is interesting to note that the deflections, moments, etc., ob- 
tained by unloading from the maximum load and subsequent loading 
to -80 psi (-5.52 x 10^ N/m 2 ) are virtually the same as those 
that would be obtained simply by loading monotonically to -80 psi 
(-5.52 x 10^ N/m 2 ) from the initial state. It is conjectured that 
this occurs as a combined result of assuming elastic-perfectly 
plastic behavior, neglecting the effects of geometric nonlinearity, 
and the fact that the same material properties were assumed to 
exist in reversed loading. Moreover, the values of residual stress, 
strain, and deflection obtained at the end of one full cycle are 
virtually the negatives of those values obtained by unloading to 
zero load from the maximum load. 


Simply Supported Circular Plate . To investigate the generality 
of these results, a different structure, a simply supported circular 
plate subjected to a uniform pressure applied centrally over a cir- 
cular area with radius 0.0718 of the plate radius, was cycled 
through various load ranges. Again, elastic-perfectly plastic be- 
havior was considered. The material properties assumed were E = 

10.5 x IQ 6 psi (7.24 x 10 10 N/m 2 ), v = 0.33, a Q = 4000 psi 

(2.76 x 107 N/ m 2) # The radius of the plate was 2.61 in. (6.63 x 
10*2 m ) s and the thickness was 0.2615 in. (6.64 x 10"^ m) . The 
load ranges considered were ±2000 psi (1.38 x 10 7 N/m2), ±3000 psi 
(2.07 x 107 N/m 2 ), ±3500 psi (2.41 x 10 7 N/m 2 ), and ±4000 psi (2.76 x 

10 7 N/m 2 ). The results are presented in Fig. 15. In all cases 
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except the last, the displacements, moments, etc., at the maximum 
negative load obtained by unloading from the maximum positive load 
are the same as would be obtained merely by loading monotonically 
to the maximum negative load from the virgin state. For the last 
case [±4000 psi (±2.76 x 10 7 N/m 2 ) load range] the load incre- 
ments used during reversed loading are too large from the stand- 
point of accuracy, and consequently, the plastic strains computed 
are smaller than those that actually occur [4000 psi (2.76 x 
1G 7 N/m 2 ) is near the theoretical collapse load of 4280 psi 
(2.95 x 10 7 N/ni ) for this structure]. These cases tend to cor- 
roborate the hypothesis that for elastic-ideally plastic material 
one need only consider one-half cycle of loading to obtain informa- 
tion concerning full cycle behavior when the effects of geometric 
nonlinearity are ignored. 


Clamped Circular Plate . A strain hardening problem is con- 
sidered next. A uniformly loaded clamped circular plate was 
cycled between ±560 psi (±3.86 x 10^ N/m 2 ). The same problem 
was considered for monotonic loading up to 560 psi (3.86 x 
106 N/m 2 ) in Ref. 49, and the results for this range are com- 
pared. Excellent agreement up to the maximum load was achieved 
(see Fig. 16). The discrepancies at this load may be attributed 
to the use of different plasticity theories (kinematic versus iso- 
tropic hardening) and the difficulty in reproducing the stress- 
strain data from Ref. 49. Furthermore, the load-deflection curve 
does exhibit all of the characteristics of strain hardening be- 
havior. The absolute magnitude of the center deflection at the 
maximum negative load is larger than that developed at the maxi- 
mum positive load, and the full cycle residual deflections axe 
triple those of the half cycle. 


Stiffened Shallow Spherical Shell . Results for a uniformly 
loaded shallow spherical shell with a stiffened circular hole at 
the apex are presented next. This problem demonstrates the bene- 
ficial effects of a stiffening ring on the elastic-plastic behavior 
of a shell, although for this particular problem it is seen that 
large deflection terms are also important and should be included. 
The pertinent geometric and material parameters defining the prob- 
lem are shown in Figs. 17 and 18. The ratio of the hole radius to 
shell base plane radius (b/a) is 0.1, and elastic-perf ectly 
plastic material behavior was assumed for the shell. 


Figure 17a shows the normal displacement versus the applied 
pressure at the ring hole interface and at an interior point ap- 
proximately halfway between the hole and the outer edge boundary 
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for an unstiffened hole and one with a stiff ring. As seen from 
Fig. 17a there is a substantial difference between the displace- 
ment at the hole boundary for the case of an unstiffened and for 
that of a stiffened hole. In fact, at the collapse load the dis- 
placement for the stiffened hole changes sign and is in the di- 
rection opposite to that of the applied uniform pressure. In ef- 
fect, the region in the vicinity of the hole moves as a rigid body 
as the displacements in the interior become unbounded. This is 
due to the restraining effect of the ring in preventing the bole 
circumference from contracting. Since the effect of the hole is 
localized, the displacements in the interior (Fig. 17b) for the 
stiffened and unstiffened case are indistinguishable. As indi- 
cated, sudden collapse of the shell is evidenced at qa^/Et^ - 15000. 
This occurs when the entire cross section in a substantial portion 
of the interior is plastic for both cases considered. However, 
since the ring carries a portion of the load, there is a wholly 
elastic section between the hole boundary and completely plastic 
Interior cross section at collapse. This contrasts with the un- 
stiffened case, for which the wholly plastic cross sections begin 
at the hole boundary and propagate towards the interior with in- 
creasing load. 

Figures 18a and b show the distribution of circumferential 
stress resultant at the yield load and at an intermediate load in 
the plastic range. As expected (Fig. 18a), the peak value for the 
unstiffened hole is at the hole boundary. As the region of plas- 
ticity expands, this peak value moves toward the interior and is 
located approximately at the elastic-plastic boundary. Figure 18b 
shows results at the same two loads for the stiffened hole. It can 
be seen that the stiff ring substantially reduces the stress re- 
sultant at the hole boundary. 


Analysis of Three Dimensional Bodies 

The isoparametric hexahedra family of elements presented in 
Ref. 50 represents the element library for this program. 

End-Loaded Rectangular Prism . Results for the longitudinal 
and transverse stress distributions in a rectangular prism sub- 
jected to prescribed end forces, shown in Figs. 19a and b, are 
used to demonstrate the accuracy of this element for elastic be- 
havior. The idealization used to represent a quadrant of the 
prism consists of 128 members and 225 nodes and resulted in 
560 degrees-of-freedom when the simple 8-node linear displace- 
ment field element was used throughout the idealization. The re- 
sults are compared with a three dimensional elasticity solution 
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presented in Ref. 51. As seen from these figures the correlation 
is quite good although it should be noted that the elasticity solu- 
tion corresponds to points along the x-axis with y = 0 and 
z = 0, whereas the results from the finite -element analysis are 
plotted at points corresponding to the centroid of those elements 
adjacent to the axes of symmetry. 

Uniformly Loaded Rectangular Laminate . The elastic stress 
distribution in a laminated fibrous composite is studied by means 
of the three dimensional finite-element module of the PLANS pro- 
gram. Of interest are the normal and shear stress components in 
the interlaminar region, which cannot be predicted by using classi- 
cal plate theory. The influence of interlaminar deformation on the 
delamination and general failure of laminates is currently under 
investigation. 

The stress distribution in each layer of a boron-aluminum 
laminate with a symmetric [ 90/90/0/0] g layup is shown in Fig. 20. 
The interlaminar normal (a z ) and shear (Ty Z ) stresses are seen 
to be significant only in the vicinity of the free edges. This 
effect was noted in several previous studies, and it has been shown 
that the interlaminar stresses decay beyond a distance that is of 
the order of the thickness of the laminate. 

For the layup and loading considered, the interlaminar normal 
stress is a favorable compressive stress tending to retard de lamina- 
tion. A reversal in applied stress or layup (e.g., [ 0/0/90/90 ] s ) 

would result in an unfavorable tensile normal stress. 

Infinite Plate Containing a Circular Hole . The classical 
problem of an infinite plate containing a circular hole is chosen 
to demonstrate the use of the three dimensional element for a 
problem involving material nonlinearity. The thickness- to-hole 
diameter ratio was taken as 0.75 so that a plane stress analysis 
would not be valid [52]. The finite-element model consisted of 
four layers of eight-node hexahedra elements to represent the 
half-thickness , with each layer consisting of sixty elements. 

Results of this analysis are shown in Fig. 21 for the varia- 
tion of circumferential stress oq , along the axis of symmetry 
perpendicular to the load at the maximum elastic load. The elas- 
tic distribution along the middle and upper surfaces compares 
quite favorably with published results [52]. The values of stress, 
obtained from the finite-element analysis, at the edge of the hole 
are obtained by extrapolation, since centroidal stresses are used 
exclusively In this analysis. 


33 



The dashed curve in Fig. 21 represents , the distribution of 
oq at almost twice the maximum elastic load, for elastic-perfectly 
plastic behavior, and reveals a substantial reduction in peaking of 
the stress field in the vicinity of the hole. Although the plate 
yields initially at the middle surface, the elastic-plastic boundary 
quickly propagates through the thickness and thereafter follows a 
pattern that is similar to the behavior associated with a plane 
stress state. 

Plane Stress Analysis of Laminated Composites 

A finite -element to model the interlaminar behavior of multi- 
layered fibrous composites under plane stress loadings has previously 
been presented and discussed [53,54]. The idealized model separates 
the membrane and interlaminar properties of a laminated composite by 
using alternating orthotropic fiber-bearing segments, which are as- 
sumed to remain elastic throughout the entire load history, and iso- 
tropic, elastic-plastic shear segments. The constant stress triangle 
is used to represent the orthotropic segments that carry in-plane 
stresses only; the shear segments carry only interlaminar shear 
stresses, and are in a state of pure shear. 

Results for several sample problems are presented and discussed 
in detail in Ref. 54. One such problem involves a flat panel with 
a circular cutout, loaded along two opposite edges. The dimensions 
of the panel were chosen such that the stress field around the cut- 
out is not influenced by the external boundaries. A typical idealiza- 
tion is shown in Fig. 22, which contains 223 members and 138 nodes. 

A fine network of elements are used around the edge of the hole 
since the interlaminar region in which the stress fields diverge 
from the classical plate solution comprises a narrow region along 
the stress free edge of the laminate. The interlaminar shear stress 
distribution around the edge of the hole is shown in Fig. 23 at the 
maximum elastic load (cr 0 /Tyi e ]_cl = 2.95), maximum and minimum loads 
in the plastic range (t 0 /ry:L e ;LcL = ±7.37), and at a load corre- 
sponding to complete unloading from the minimum load. As indicated 
in the figure, high residual shear stresses occur for 6 < 45 °, and 
relatively low residual shear stresses occur for 6 > 45 °. 


Geometric Nonlinearity 


Circular Plates 

The capability of treating geometric nonlinearity and combined 
geometric and material nonlinearity has been developed under the 
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present contract for axisymmetric shell structures. The analysis 
uses the shell element reported in Ref. 42, and a presentation of 
the application of the element for combined nonlinear analysis is 
presented in Ref. 23. The accuracy of the procedure for geometric 
nonlinearity in the case of purely elastic behavior is demonstrated 
in Figs. 24 and 25 where a comparison of results obtained from the 
present analysis is made with those obtained in Ref. 55 for a 
clamped, uniformly loaded, elastic circular plate. Poisson’s ratio 
was chosen to be 0.3. Figure 24 is a plot of central deflection 
versus load, and Fig. 25 is a plot of bending and membrane stresses 
at the center and edge versus deflection. For the increment size 
chosen, excellent agreement was obtained between the solution pre- 
sented in Ref. 55 and our numerical results. 

In the present investigation, results for this problem were 
obtained by using both the "tangent modulus" method and the "effec- 
tive" load method for the same increment size. For the latter case, 
the stiffness matrix was re-formed every five increments. No equi- 
librium correction term was included for either method for this 
problem. The deflections and bending stresses in both cases were 
identical, while slightly smaller membrane stresses were predicted 
by the effective load method. Of most significance was the re- 
duction in CPU time from 386.28 seconds for the tangent modulus 
method to 202.08 seconds for the effective load method, an ap- 
proximately 47 percent time savings at no appreciable loss in 
accuracy. Similar time savings of from 40 to 50 percent were 
noted for other problems. 


Spherical Cap 

Figure 26 illustrates the need for the equilibrium correction 
term in problems involving a high degree of nonlinearity. An exact 
load-deflection curve obtained from Ref. 20, based upon results 
presented in Ref. 56, is shown for an elastic, clamped spherical 
cap loaded by a central concentrated load. Also shown are results 
obtained from the current analysis using a straight incremental ap- 
proach with 1/8 lb (0.56 N) Increments and an incremental-plus- 
equilibrium correction solution using 1 lb (4.45 N) Increments. 
The results obtained in the current analysis are virtually identical 
(for both incremental and incremental-with-equilibrium correction 
solutions) to the numerical results given in Ref. 20. 
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Combined Material and Geometric Nonlinearity 
Circular Plate 

Figure 27a shows a load versus central deflection plot for a 
centrally loaded, simply-supported plate with a diameter-to- 
thickness ratio of 40.6. The numerical results are compared with 
test data obtained from the experimental program described in 
Ref. 23. Shown are the linear elastic, nonlinear elastic, elastic- 
plastic and combined nonlinear predictions using the tangent modulus 
approach with the incremental and incremental-with-equilibrium cor- 
rection solution procedures. Although for these combined problems 
the equilibrium correction affords a considerable improvement over 
the incremental approach, without equilibrium correction, a more 
extensive iteration scheme is probably needed to close the theoretical- 
experimental gap. In Figs. 27b and 27c the radial distribution of 
circumferential strain at the lower and upper surfaces, respectively, 
for this plate is illustrated for several load levels and compared 
with theory. Despite the only fair-to-good correlation of the dis- 
placement data at high loads, excellent correlation with experiment 
for the strains is noted, except, as might be anticipated, directly 
under the loading rod for higher loads. The discrepancy at this 
point might be the result of local shear and penetration effects. 


Spherical Cap 

Wilkinson and Fulton [57] have presented results for the 
elasto-plastic buckling of uniformly loaded shallow spherical caps 
with both simple and clamped support at the edges. Several compari- 
son test cases were chosen to verify their results and determine 
the present program* s ability to predict buckling loads of such 
structures. The cases run were for a = 0.1, (3 = 0.002, and 
A = 4 and 5.5 for the clamped cap, and A = 4 for the simply- 
supported cap. Here a is the ratio of tangent modulus to Young’s 
modulus, P is the ratio of yield stress to Yyung* s i modulus , and A 
Is the geometric shell parameter 2 [ 3 ( 1 - v2)]t(H/h)2, with v be- 
ing Poisson's ratio, H the maximum shell rise, and h the shell 
thickness. As can be seen from Fig. 28, excellent agreement for 
these cases was obtained. The buckling pressures are 3 percent 
higher than those predicted by the analysis of Ref. 57. The results 
from the previous analysis are probably more accurate than those ob- 
tained here, since the present analysis makes no attempt to refine 
the load increment size in the vicinity of the critical load. The 
tangent modulus method was used to account for the effects of geo- 
metric nonlinearity. 
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Cyclic Loading of Simply-Supported Circular Plate 

Results in the form of load versus central deflection of a 
simply-supported , centrally loaded mild steel circular plate are 
shown in Fig. 29. These results involve a history of loading to 
a maximum load in the plastic range and then the removal of the 
load. A comparison with the experimental data presented in Ref. 58 
indicates that the finite-element results predict larger displace- 
ments than those obtained experimentally. This may be partially 
explained by the fact that no information (except the yield stress) 
was available in Ref. 58 concerning the strain hardening proper- 
ties of the material used in the experiment. The finite-element 
analysis was performed by assuming elastic-ideally plastic be- 
havior, which is a good representation of the stress-strain be- 
havior for mild steel for strains less than 2 percent. The 
larger displacement prediction from analysis is consistent with 
this assumption. When the strains in the plate become larger than 
2 percent, mild steel experiences strain hardening. Indeed at 
loads above 15,000 lb (66.73 x 1C' 5 N) the theoretically pre- 
dicted strains exceed 2 percent in a considerable region of the 
plate, and divergence of the results occurs. As a consequence of 
the overprediction of the maximum displacement, the residual dis- 
placement predicted by the analysis is considerably greater than 
that experimentally observed. However, the general shape of the 
load-deflection curve upon unloading parallels the experimental 
curve . 


Ring-Stiffened Cylinder 

The circular cylinder is a structural component that is cur- 
rently under consideration as an energy-absorbing device to improve 
the crashworthiness of small land and air vehicles. A reasonably 
accurate prediction of the load-deformation behavior of this de- 
vice is required to evaluate its efficiency as compared to other 
devices of a similar kind. Toward this goal the prototype program 
was used to describe the elastic-plastic, large deflection behavior 
of a simply-supported ring-stiffened circular cylinder shown in 
Fig. 30. 

Results for the axial displacement versus axial stress resultant 
(Fig. 30a) and the radial displacement (Fig. 30b) are shown for elas- 
tic and elastic-ideally-plastic , large deflection behavior. The 
rings (EA *= 10& lbs = 4.45 x 10^ N) remain elastic throughout the 
entire load history, but the influence of material nonlinarity in the 
cylinder significantly influences its behavior and decreases Its 
buckling load. From the results for the radial displacements, it is 
apparent that the rings have a negligible influence on the behavior 
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of this cylinder. To verify this conclusion, the problem was re- 
analyzed (elastic and elastic-plastic) for the unstiffened cylinder, 
and it was found that the removal of the rings had virtually no ef- 
fect on the buckling load or mode shape. 

Unstiffened Circular Cylinders 

Results for the end-load versus end-deflection for a clamped 
circular cylinder fabricated from 6061-0 aluminum alloy are pre- 
sented in Fig. 31a. They indicate a buckling load of about 
1600 Ib/in. (280 x 10 2 N/m) or 39,800 pounds (177 x 10 2 N) 
total load for this structure. This compares quite favorably 
with experimentally determined buckling loads of 35,650 pounds 
(159 x 10 2 N) (two specimens) and 41,250 pounds (183 x IQ 2 N) 

(one specimen) obtained at NASA/Langley. In the course of the 
analysis , it was found that the buckling load is extremely de- 
pendent on an accurate determination of the material properties. 

A Ramberg -Osgood three parameter representation of the stress- 
strain curve for 6061-0 aluminum alloy was used for the analysis, 
i . e , , 
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of 4300 psi (2.97 x 10 7 N/m 2 ) was used. These properties were 
derived from data obtained from the Materials Selector Handbook 
and verified by tests performed in the Grumman Research Department 
Applied Mechanics Laboratory. Pertinent geometric data for the 
clamped circular cylinder are given in Figs. 31a and b. Calcula- 
tions were then performed, using the same configuration and mate- 
rial properties, but with simple supports at the ends to determine 
the effect of edge restraint on the buckling load of this cylinder. 
As anticipated, the reduction in buckling load resulting from the 
change in edge restraint was extremely small and virtually negli- 
gible (see Fig, 31a). 


Figure 31b shows the radial deflection profile versus axial 
length for the clamped and simply- supported cylinders at end-loads 
of 1050 Ib/in. (184 x 10 2 N/m) and 1575 Ib/in. (276 x 10 2 N/m), 
which is immediately before buckling. The initial buckling pro- 
files are quite evident, with the wavelength of the simply-supported 
cylinder being smaller than that of the clamped cylinder, and the 
maximum displacement amplitude of the simply-supported shell ap- 
proximately 1 \ times that of the clamped cylinder at this load, 
indicative of the slightly lower critical load for the simply- 
supported cylinder . 



Clamped. Truncated Circular Cone 

A clamped truncated circular cone was then analyzed in an ini- 
tial effort to determine what qualitative and quantitative bene- 
ficial effects tapering the sides might introduce into the energy 
absorption capacity of a vehicle. The first cone chosen had the 
dimensions of one specimen of the test program at Langley. It, 
too, was fabricated from 6061-0 aluminum alloy and its dimensions 
are presented in Figs. 32a and b. The cone angle, a, is 17.09 de- 
grees. The critical load for this structure was calculated to be 
42,575 pounds (189 x IQ 3 N) , corresponding to an axial stress re- 
sultant applied at the larger end of 670 Ib/in. (117 x 1Q3 N/m) . 

This represents an increase of almost 7 percent over the cal- 
culated critical load for the clamped circular cylinder previously 
analyzed. The cone has the same thickness and axial length as the 
cylinder and its smaller diameter is equal to the diameter of the 
cylinder. Figure 32a shows an applied load versus end-deflection 
curve for this structure, and Fig. 32b presents a normal displace- 
ment profile versus axial distance for loads of 420 lb/in. (74 x 
10 3 N/m), and 660 Ib/in. (116 x 10 3 N/m) (i.e., immediately pre- 
ceding buckling) . It can be seen that the structure buckles at the 
smaller diameter end. 

Several additional analytical cases have been considered because 
the experimental buckling loads obtained for the cone with a = 17° 
were considerably lower than those obtained analytically assuming 
clamped edges. It was observed during the experiment that, at the 
small diameter end, the edges "curled" as the load increased,, in- 
dicating that clamped edge conditions are not justifiable. To deter- 
mine the extent of the effect of edge restraints, the same cone was 
analytically buckled with the small diameter end simply supported and 
free (axially restrained with no radial force or bending moment al- 
lowed) . The critical load for the simply supported case was slightly 
lower than for the clamped case and was equal to 38,800 lbs (172.7 x 
10 3 N) . For the left end free, the buckling lead was drastically re- 
duced from 42,575 lbs (189 x IQ 3 N) , for the clamped case to 
13,000 lbs (57.9 x 1G 3 N) , a 69.4 percent change. (See Fig. 32a 
for load-end deflection curves for the simply supported and free end 
cases.) The experimentally observed actual load of 25,000 lbs 
(111 x XG 3 N) lies between the two extremes of edge restraint repre- 
sented by simple supports and a free edge, thus indicating the im- 
portance of accurately reproducing the restraint conditions for the 
plastic buckling analysis of conical shells. 
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4. PROGRAM DEFINITION 


As mentioned in Section 1, the computer program PLANS for the 
PLastic ANalysIs of Structures is organized as shown in Fig. 2b. 

The basic philosophy of program organization is to separate classes 
of analysis into individual groups, with each group defined by the 
physical problem to be solved. For example, as indicated in the 
figure, groups can be defined for membrane stressed structures 
(PLANE), combined bending and membrane stressed structures (BEND), 
bodies of revolution (REVBY) , three dimensional solids (HEX), and 
membrane stressed laminated composites (COMPEL) . On the basis of 
this concept, each group is in itself an independent finite-element 
computer program, with an associated element library, that can be 
individually loaded and used to solve the problem class of interest. 
In this manner any simplification or specialization germane to the 
individual analysis can be incorporated. Each group can contain 
subprograms common to all other program groups as well as subpro- 
grams distinct to its group . Linkage to any individual program 
group is made through an executive program indicated in the figure. 
As shown, linkage to any program group can be made only through 
this executive routine with no horizontal linkage between group 
programs . 

A brief discussion of the significant features of several 
components of the program follows. A more detailed discussion 
is presented in a separate report. 


Executive Program 

The executive program performs the function of loading the 
individual analysis programs. Data set management and core al- 
location are handled by each of these programs. The system allows 
a user to store many program decks on a magnetic tape or direct 
access volume, and to use these programs conveniently. During a 
single computer run, the user may call foir a n Select-Load~Go" pro- 
cedure to execute routines in a previously created object file, or 
he may issue a series of easily-learned commands to change the old 
source file, thereby generating new source, index, and object files 
incorporating the desired changes. These new files can then be 
made available for a "Select-Load -Go" as part of the same run. 

The executive system comprises three separate collections of 
subprograms to accomplish this. These are: 

® A source update program to allow a user to 
selectively edit and compile a new source 
program. 



® An object update program that maintains up- 
dated files of compiled source programs. 

® A file select program that selectively loads 
programs from the object file for execution. 

While the file select program can accept commands to load in- 
dividual programs, it can also load groups of programs bjr simply 
supplying a group name that has been previously defined. In addi- 
tion, each group definition may contain other group names in Its 
definition. This feature is particularly meaningful for our pur- 
poses since each program in Fig. 2b naturally defines a group, and 
programs common to all groups (such as the solution package) are 
naturally subgroups. Thus, as shown in Fig. 2b, groups can be set 
up corresponding to PLANE, BEND, HEX, etc., and subgroups can be 
set up defining, for example, the solution package. 

Thus, the capability to maintain source and object files and 
selectively load programs by groups is the most significant feature 
of the executive program. In addition, this system will have the 
ability to update files by rereading complete groups and is able to 
define new groups as additional programs are added to the PLANS 
system. 


Basic Flow of Analysis Programs 

Figure 33a is a gross schematic representation of the flow of 
each of the analysis programs. As shown in the figure, each module 
has three components: a main calling routine (MAIN); and elastic 

analysis group (ELAS) ; and the plastic analysis group (PLAS) . A 
brief description of these components follows. 


Main Program 

Each individual computer program defining a group is controlled 
by a main calling program. This program sets up core allocation for 
principal arrays and the data set specifications for auxiliary 
storage. It then transfers control to subroutine ELAS. ELAS is a 
finite -element elastic analysis program that performs the elastic 
analysis and calculates the initial yield load. Control Is then 
transferred back to the main program, which calls Subroutine PLAS 
only if requested. PLAS manages the plastic analysis and main- 
tains control of the analysis until the complete plastic analysis 
Is performed. At this point control Is transferred to the main 
program either to terminate the job or call ELAS again to work 
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another problem. By so organizing the program it should be possi- 
ble with a minimum amount of changes to use PLAS with other avail- 
able finite-element elastic analysis programs. 


ELAS 


Figure 33b shows a block diagram of the computational flow of 
ELAS. This program is a special purpose finite “element program for 
the elastic analysis of structures. Accordingly, its first major 
task is to read all inputs. The input is read in functional groups 

as follows: 

• Problem title. 

® Nodal coordinates and control variables. 

• Member topology describing connectivity. 

® Nodal boundary conditions. Single and multi- 

point constraints on the displacements are 
specified in addition to fixed or free condi- 
tions . 

® Load vector — includes the consistent member 
load vector for a distributed load. 

® Material properties — tables of elastic and 

plastic material properties and member 
geometric properties are set up along with 
applicable members. Complete generality has 
been maintained so that orthotropic material 
behavior can be specified. 

Other variables controlling output are also specified. The input 
scheme has been written so as to minimize the amount of card 
handling for any problem. 

The next step is to form all element stiffness, stress, and 
initial strain stiffness matrices. This routine also places the 
elements of the stiffness matrix along with their position in 
the total stiffness matrix on an auxiliary storage device. 

The total stiffness matrix and load vector are assembled by 
sequentially reading this information, stacking that portion of 
the stiffness matrix that fits in core, and then reading it onto 
an auxiliary storage device. This process is repeated until the 
entire load vector and stiffness matrix have been formed. 
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At this point it is appropriate to mention that our basic de- 
sign philosophy is to make alterations, perform certain matrix 
multiplications, account for boundary conditions on the element 
level, and then to assemble these quantities into total arrays. 

In this manner we need not make use of a matrix interpretive sys- 
tem to manipulate large matrices. 

With the total stiffness matrix and load vector assembled, the 
next step is to solve for displacements. Two solution packages 
are currently being used for this purpose, both based on the Cholesky 
factorization method . Their differences are based on the system 
used to store the stiffness matrix. PIRATE is based on partition- 
ing the structure so that the stiffness matrix is explicitly in 
block tridiagonal form. The second routine, PODSYM, is based ex- 
plicitly on the banded nature of the stiffness matrix, so that only 
elements within the lower triangle that lie between the semiband- 
width need be assembled. In addition, this routine "packs" the 
rows of the stiffness matrix before writing it on auxiliary storage 
by suppressing all consecutive zeros. 

The unit load stresses are calculated next from displacements. 
These stresses are used to calculate the initial yield load, and 
this yield load is used to scale the displacements and to calculate 
yield load stresses and strains. Control is then returned to the 
main program. 

Subroutine PLAS 

Figure 33c shows a block diagram of the computational flow of 
Subroutine PLAS. This program supervises the entire plastic analy- 
sis after initial processing has been carried out by ELAS. The 
principal information that must be communicated to PLAS is: 

® Factored stiffness matrix and unit load vector 

® Element stress and initial strain matrices 

® Initial yield load 

® Plastic material properties. 

With this information, the routine increments the load and 

checks the yield condition for each member. Members that were pre- 
viously plastic are also checked against an unloading criterion. 

The effective plastic load vector is then formed as the product 
of the member’s initial strain matrices and the plastic strains and 
added to the vector of applied load. This combined load vector Is 
then used as the right hand side of the matrix equation to calculate 
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nodal displacements. It should be noted that, at this point, the 
stiffness matrix has been factored so that the solution involves 
only the multiplication necessary to perform the forward and back- 
ward solution. 

In the next step, a subroutine is called that implements the 
plastic constitutive relations. This segment of code, calculates 
increments of total strain, stress, and plastic strain (shift in 
the yield surface when kinematic hardening is used) for each member 
from the increments of displacement, and sums the results with 
total quantities calculated previously. 

After control is returned from the above subroutine a test is 
made to see if the maximum specified load has been reached. If not, 
the load is incremented and the steps discussed above are repeated. 

When the maximum load is reached, the final effective plastic 
load vector is formed and displacements calculated with the applied 
load vector set equal to zero. Stresses and strains are calculated 
on the basis of these displacements. If subsequent yielding occurs 
in the reversed cycle at a load that is opposite in sign to the load 
previously completed, then these stresses and strains represent 
residual quantities. This is checked next by computing a new yield 
load. This calculation involves the solution of a quadratic equa- 
tion for each member. One root of this solution represents the 
previously reached maximum load and the other the new yield load in 
the reversed direction. This load level is used as a starting 
point for the next cycle of loading. 

At this point, a check is made to see if another cycle of load- 
ing is desired. If not, control is returned to the main program. 

If an additional cycle of loading is specified, new plastic mate- 
rial parameters can be read as input and inserted into the material 
property tables. This is followed by the formation of the effec- 
tive plastic load vector. This load vector is then added to the 
load vector, and the sum is used to solve for the new critical load 
displacements and the corresponding stresses and strains. Transfer 
is then made to the beginning of PLAS in order to increment the 
load and proceed as described above. 


Representative Time Requirement s 

The cost involved in obtaining an elastic solution to a par- 
ticular problem depends upon a multitude of parameters. These in- 
clude the core requirements, number of operating system calls, the 
total and central processing time, etc. An attempt to uniquely de- 
fine the cost requirements for a problem of a given size, run under 
any one of the current operating systems, is unfeasible. 



Although a realistic accounting procedure would consider some 
or all of the variables mentioned, it nevertheless remains true that 
the central processing time (CPU) is a significant factor in deter- 
mining the total cost. On this basis we have indicated (Fig. 34) 
representative elastic solution time requirements versus problem 
size for the PLANE module. The problem size is represented in the 
figure by the number of degrees-of-freedom and the semi-bandwidth. 
The two curves in Fig, 34 were not obtained independently. That 
is, as illustrated in the figure, a problem involving 700 degrees- 
of-freedom requires 13 seconds of CPU if the semi -bandwidth is 90. 

The time requirements for the solution to an elastic -plastic 
problem is measured as the summation of the time required for each 
increment of load. As discussed in Section 2, the displacement 
method- predictor procedure (used by all the modules of PLANS) does 
not require the reformulation and successive decomposition of the 
elastic stiffness matrix for each step beyond the initial maximum 
elastic load value. Thus, the time required for each incremental 
load step is considerabljr less than the elastic solution require- 
ment. In general, it has been found that one-tenth of the elastic 
solution time represents a fairly reliable estimate for determin- 
ing the CPU time for each increment in the plastic range. 

A schematic curve of total system time versus problem size 
(degree-of-freedom x semi-bandwidth) is presented in Ref. 59, in 
addition to an excellent discussion on the cost of computing for 

a nonlinear analysis system. 


5. CONCLUDING REMARKS 


The material presented in this report, although primarily con- 
cerned with the authors’ investigations, is representative of the 
current state of the development of nonlinear analysis techniques 
within the framework of the finite-element method. In addition to 
this report, recent surveys of the literature reveal that these 
methods have been developed to the point where a broad spectrum of 
structures of considerable complexity can be analyzed through the 
plastic range. Particular emphasis has been placed here on cyclic 
loading conditions involving reversed loading into the plastic 
range. Information concerning the behavior of structures subjected 
to realistic load spectra of this type should be of great value in 
treating the problem of low-cycle fatigue. 
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Although deficiencies still remain in the development of ac- 
curate constitutive relations describing elastic-plastic large 
deformation behavior, we do not believe that this should inhibit 
the development of analytical methods and their associated com- 
puter programs. Whether any one of the currently available plas- 
ticity theories represents an "exact" model for material nonlin- 
earity is of less significance than the fact that the finite- 
element method allows us to perform numerical experiments to aid 
in the development of theories that more realistically describe 
complex material behavior. The programs also provide a guide for 
the experimentalist so that he can determine the type of experi- 
ments to perform, the properties to measure, and the data to 
monitor so that more accurate theories may be formulated. On this 
basis it is anticipated that further inroads can be made into those 
areas of investigation that, in the authors 5 opinion, require fur- 
ther insight. These problem areas include: anisotropic plasticity, 

temperature dependent plasticity, and dynamic loading effects. 

During t he course of the present study a substantial effort 
was devoted toward the development of a user-oriented computer 
code, briefly described in this section, that implements the plas- 
tic analysis techniques favored by the authors. To translate 
these theoretical approaches into a user-oriented, practical 
computing tool, one must recognize that those problems that exist 
in the development of such a tool for an elastic analysis are 
magnified for a nonlinear analysis. The constraints imposed by 
the range of validity of any particular theoretical analysis must 
be understood by the user and the data obtained from it carefully 
interpreted to ensure its proper application. 'For these reasons 
it would be potentially dangerous to consider such a tool as a 
"black box." Thus, a desirable feature of such a program, i.e., 
its ease of application, perhaps represents its greatest liability. 

Recently, considerable attention [18] has been turned toward 
considering problems of efficiency, accuracy, and the resultant 
economic feasibility of large scale use of numerical nonlinear 
analytic systems. On the basis of the improvements in solution 
procedures, equation solving algorithms and computer hardware, 
one can now envision the time when it will become unnecessary to 
guess at or design around a problem involving complex nonlineari- 
ties . In considering the cost of using such a tool, however, one 
must always weigh the cost of the analysis versus the cost and 
risk involved in the possibility of designing a system on the basis 
of partly al information or instituting an experimental program to 
gain the same information. 
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APPENDIX A 


CATALOG OF FINITE-ELEMENTS 


Plane Stress (Membrane States Only) 


1. Constant Strain Triangle (CST) (Fig. 35) — This well" 
known plane stress membrane element was used successfully for the 
idealization of structures considered under Contract NAS 1-5040 . 

Its derivation is based upon the assumption of a linear distribu- 
tion for the in-plane displacements u and v, and consequent!}?, 
leads to a constant strain state within the element. Each vertex 
is allowed two degrees of freedom (the in-plane displacements u 
and v) for a total of six degrees of freedom for the element. 
Consistent with the total strain distribution, the initial strains 
(plastic strains) are assumed to be constant within each element. 
Stiffness and initial strain matrices have been developed and suc- 
cessfully, used in Contracts NAS 1~5040 and NAS 1-7315, Refs, 25 and 
16. 

2. Linear Strain Triangle (LST) (Fig. 36) — In regions of 
high strain gradient, the CST triangle is not sufficiently accurate 
to be used in a plasticity analysis unless a very fine grid is em- 
ployed. The linear strain triangle (LST) remedies this shortcoming 
of the CST element. The assumption of a quadratic distribution for 
the in-plane displacements allows for a linear strain variation 
within the triangle. Two degrees of freedom at each node (u,v) 
for each of the six nodes (three vertex and three midside nodes) 
give this element a total of 12 degrees of freedom. The initial 
strains are assumed to be constant within each element and are 
evaluated at the centroid. Both stiffness and initial strain 
matrices have been developed and successfully used in Contract 

NAS 1-7315. 

3. Hybrid Triangles — In transition regions, i.e., regions 
in which stresses and strains change from rapidly varying to slowly 
varying, it becomes convenient and efficient to switch from linear 
strain triangles to constant strain triangles. This is accomplished 
by using four and five node triangles to maintain compatibility 
with both the CST and LST elements. These elements together with 
the CST and LST elements were originally used in Ref. 6, and are 
referred to as the TRIM 3 through TRIM 6 family. For these mixed 
formulation hybrid elements, the displacements along edges may vary 
quadratically or linearly, depending upon whether an LST or CST tri- 
angle is contiguous to the respective sides. Again, the plastic 
strain distribution is assumed constant within each element. 
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4. Warped Quadrilateral Shear Panel (Fig. 37 ) — For many air- 
craft structures the stiffness properties may be adequately repre- 
sented by an assemblage of beams, bars, and shear panels. The 
shear panels, as their name implies, provide resistance only to 
shearing forces. The normal and bending stresses are carried by the 
bars and beams. The shear panel used in the analysis is based upon 
the original formulation and approximation proposed by Garvey 

(Ref. 37) . It is quadrilateral in planform and all four nodes need 
not lie in a plane. The out of plane "twisting" is taken up by 
"kick" forces fj_ as shown in Fig. 37. The plastic behavior of 
this element is assumed to be uniaxial, i.e., the single stress or 
strain component at the centroid of the element is used to compute 
the plastic strain at that point. 

5. Stringers (Fig. 38) — For many aircraft structures, e.g., 
fuselages, wings, etc., local stiffening is required to provide 
adequate stability in compression. Special one dimensional finite- 
elements are required to represent the stringers used for this pur- 
pose. Uniform cross section stringer elements have been developed 
using both constant and linearly varying strain assumptions, so 
that comnatibility with both the GST and LST elements can be main- 
tained. For the constant strain stringer, a linear axial displace- 
ment and a constant initial (plastic) strain distribution within 
the element are assumed. The linear strain stringer stiffness 
matrix is based upon a quadratic displacement assumption and the 
initial (plastic) strains are constant within the element 

Bending 


1. Triangular Element (Bell) (Fig. 39) — While a rectangular 
element is excellent for rectangular regions, a triangular element 
Is necessary to treat plates with arbitrary boundaries, including 
curved boundaries . The fully conforming triangular bending element 
developed by Bell (Ref. 39) and several other authors (Refs. 60 and 61) 
has proved to be extremely effective for the analysis of bending 
of plates with curved boundaries. This element was used effectively 
in Ref. 16 to treat triangular, circular, annular, and rectangular 
plates. The lateral displacement w is assumed to be a complete 
quint ic polynomial in x and y, the in-plane coordinates. Six 
nodal degress of freedom (w, w, x ^ w,y w, xx ^ w, X y^ w,yy) at each 
vertex node and three normal slopes at the midpoints of the sides 
(w, n ) are initially prescribed. The latter three degrees of free- 
dom are then eliminated by imposing a cubic variation of the normal 
slope along the edges, resulting in an 18 degree of freedom element. 
The plastic strain distribution is allowed to vary arbitrarily in 
the plane of the element. Plastic strains are evaluated at Gauss 
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points and the plastic load vector is formed using two dimensional 
Gauss -Legendre integration in the plane of the element. The plas- 
tic strain variation through the thickness at each Gauss point is 
monitored at a given number of points through the thickness (layered 
approach) . The strains are then numerically integrated through the 
thickness for use with the initial strain stiffness matrix to form, 
the plastic load vector. For problems involving combined bending 
and stretching the in-plane displacements are assumed to vary lin- 
early, as in the constant strain membrane triangle or cubically 
(see Ref. 10) as in Felippa’s formulation. Stiffness, initial 
stress stiffness, and initial strain matrices have been derived 
and successfully employed in the work of Ref. 16. 

2. Beam Element (Fig. 40) —-The effects of local stiffening 
may be represented by a beam element that assumes the gross geo- 
metric and elastic properties of the actual stiffener. Such a beam 
element with a rectangular cross section has been developed and 
used extensively in Refs. 13 and 62 for bending, combined bending 
and stretching, geometric nonlinearity, and integrally stiffened 
plate problems. Bending, membrane, and torsional rigidity are all 
accounted for in the stiffness matrix. The lateral deflections are 
assumed to be cubic in the axial coordinate, while the axial dis- 
placement is linear. Torsional rigidity is included through the 
assumption of a linear variation of the angle of twist. 

In addition to this standard beam element, another element has 
been developed for use as a stiffener in conjunction with the tri- 
angular plate element. This element, depicted in Fig. 40, employs 
a cubic displacement assumption for both lateral displacements and 
the axial displacement. The angle of twist is assumed to vary lin- 
early from node to node. 

Plasticity is accounted for by assuming a linear variation 
from node to node of plastic strains along the beam axis and allow- 
ing an arbitrary variation of plastic strains through the cross- 
sectional area. The plastic strain distribution is integrated 
across the area using Gauss-Legendre integration and along the 
length in closed form. 

Axisymmetric Bodies of Revolution 


1. Axisymmetric Solid Element (Fig. 41) — - The analysis of 
thick solids of revolution subjected to arbitrary mechanical loads 
is of basic interest in the aerospace industry. A triangular ele- 
ment of revolution used to analyze these problems was developed by 
Wilson (Ref. 41) and others. This is essentially a two dimensional 



element with the formulation modified to account for hoop stresses 
and strains (due to the curvature of the structure) . Its derive*" 
tion is based upon a linear assumption for the radial and axial 
displacement components. A modification to this element to include 
a torsional degree-of -freedom is incorporated for the element in- 
cluded in this module. Thus, the element has three degrees-of- 
freedom at each nodal circle. Because of its triangular shape, 
it is capable of representing arbitrary bodies of revolution. The 
plastic strain distribution to be used in the formation of the ini- 
tial strain matrix is assumed constant within each element. Stiff- 
ness and initial strain matrices have been developed and incorpo- 
rated into the program. It should be mentioned that good elastic 
results have been obtained by using this element (Ref. 63), and 
Marcal (Ref. 64) and Armen, Levine, and Pifko (Ref. 65) have ob- 
tained good results in elastic-plastic analyses. 

2. Axisymmetric Thin Shell Element (Fig. 42) — For shells 
whose thickness- to-mean-radius ratio is less than 1/10 (in many 
cases this ratio may be as high as 1/3) and where local effects 
are negligible (e.g., junctures, branches, welds), the use of a 
three dimensional element becomes wasteful. For such structures, 
shell theory, which is a two dimensional approximation to three 
dimensional elasticity theory, gives an accurate description of 
the state of stress and deformation. Moreover, since thickness 
effects are included in the kinematic model, only the circumfer- 
ential and meridional variation of pertinent variables must be con- 
sidered. This results in a considerable saving with regard to 
computer storage, since a smaller number of elements is required 
to describe the structure. 

The shell of revolution is a basic component of many aero- 
space structures. Cylinders, spheres, cones, and various combina- 
tions thereof are widely used. For these reasons, an axisymmetric 
thin shell element capable of treating arbitrary shells of revolu- 
tion under axisymmetric loading is included in the plasticity pro- 
gram . 


The shell element chosen for use here is based upon the ele- 
ment developed by Levine and Armen in Ref. 42. The derivation of 
the stiffness properties of this element is based upon a shell 
theory proposed by Sanders (Ref. 66) that uses Kirchoff’s hypothe- 
sis and Love's first approximation. The meridional displacement 
is assumed to be a cubic function of the meridional coordinate £ 
(see Fig. 42) and the normal displacement is also chosen to be 
cubic in i. The actual meridional shape of the shell is repre- 
sented. by a cubic variation within the element. Plasticity has 
been accounted for by dividing the shell into layers through the 
thickness. For the current analysis, the plastic strains are as- 
sumed to vary linearly along the meridian in the derivation of the 
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initial strain matrix. Initial strain, initial stress, and conven- 
tional stiffness matrices have been developed and successfully 
employed in combined problems of material and geometric nonlinear- 
ity. 


3. Axisymmetric Thin Ring Element — A thin ring element to 
account for the effects of local stiffening on the behavior of 
axisymmetric structures has been developed. The ring may be of 
arbitrary cross section and eccentrically attached to the struc- 
ture. Torsion, extension, and bending are accounted for. The ef- 
fects of prestress, warping, and shear deformation have been 
neglected, and the shear center and centroid must coincide. The 
element is basically similar to that presented by Cohen (Ref. 43) . 
Plastic effects are included by numerically integrating the plas- 
tic strains over the cross-sectional area of the ring for use with 
the effective plastic load vector. 


Three Dimensional Analysis 

Isoparametric Rexaheuron Element Family Ranging from 8 Nodes 
(at the Corners) to 20 (8 Corner and 12 Mid-edge) Nodes (Fig. 43 
and Ref. 50) — For the 8 node element the displacement function 
is represented by a set of linear interpolating functions. Qua- 
dratic interpolation functions are used for the 20 node element. 

In local coordinates the element is a cube; in global coordinates, 
however, the position of all nodes is arbitrary. Each edge of the 
cube is either a quadratic curve in space or a straight edge, de- 
pending upon whether the edge is defined by the location of three 
or two node points . Thus, curved boundaries are represented by 
quadratic interpolating polynomials. 

A constant plastic strain within each element has been in- 
corporated. 


Laminated Composite Plate Analysis — 


® Constant stress element , 

• Pure shear element } Com P° site Element (Fig. 44) 

Interlaminar shear deformation, the mechanism by which load 
is transferred through a matrix material between two stiff laminae 
as the laminae tend to slide over each other, cannot be predicted 
by using classical plate theory. This type of deformation develops 
along the edges of a laminate and can be important with respect to 
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strength predictions of composite structures, especially for a lami- 
nate with a relatively low transverse shear strength. Analogous 
behavior is found to exist in bonded structural joints. 

The idealized model used in the analysis separates the mem- 
brane and interlaminar properties of a laminated composite by using 
alternating orthotropic fiber-bearing segments and isotropic shear 
segments, as shown in Fig. 44a, The ortho tropic segments carry in- 
plane stresses only, and may be considered to be in a state of plane 
stress; the shear segments carry only interlaminar shear stresses, 
and are in a state of pure shear. 


The composite element consistent with the idealized model is 
shown in Fig. 44b. The membrane segments are triangular orthotropic 
elements in which the total strains are assumed to be uniform. 

Here the strain-displacement relation is based on a linearly vary- 
ing displacement field. The stiffness properties of the interlami- 
nar shear segments are also based on a linear displacement field, 
so that the shear strains may be written in terms of the nodal dis- 
placements in the following manner: 


xz 


du / £-1 , £~1 , £- 

~ = u . + u . + u. 

oz \ i j k 
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where u and v are displacements in the x and y directions, 
respectively; subscripts identify elements vertices; and super- 
scripts identify element faces as shown in Fig. 44. Since the dis- 
placements vary linearly in the plane of a membrane segment, the 
interlaminar shear strain is computed on the basis of centroidal 
values of displacement, and thus the shear segment may be regarded 
as a shear-resisting medium connecting the centroids of adjacent 
membrane segments. Any number of segments may be stacked through 
the thickness to form the multilayered composite element. 


Plastic behavior of the plate is assumed to be confined to the 
softer matrix material. A uniform plastic shear strain is assumed 
to exist in the interlaminar region. A detailed discussion of this 
element and its potential usage is presented in Ref. 54. 
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APPENDIX B 


PLASTICITY RELATIONS 


Considerable efforts have been made in the experimental ex- 
ploration of the nonlinear behavior of many of the commonly en- 
countered structural materials, A motivating force underlying 
much of this effort can be attributed to basic investigations con- 
cerned with low cycle fatigue (Refs. 67 and 68). From a metallurgi- 
cal viewpoint, the mechanism associated with the behavior of duc- 
tile materials experiencing cyclic plastic deformation to fracture 
appears to be overwhelmingly complex. Because of this complexity 
no universally applicable laws governing the behavior of materials 
in the plastic range have yet been developed. If the structural 
analyst is to predict the gross deformation behavior of structures 
in the plastic range, he is required to choose from among the sev- 
eral available plasticity theories one that successfully combines 
mathematical simplicity with a reasonably faithful representation 
of some of the more obvious experimentally observed features of 
material behavior. 

We here present and briefly discuss plasticity relations for 
three particular theories associated with strain -hardening behavior. 
These are the kinematic hardening and isotropic hardening theories 
and a third theory that incorporates features of both. Isotropic 
hardening is restricted to monotonic loading conditions, whereas 
the last two are specifically developed to treat cyclic loading 
conditions involving stress reversals. In addition, consideration 
is given to those relations appropriate to ideally plastic material 
behavior . 

For isotropic, isothermal conditions, we can define a function 
of the stress components to be used in describing the limits of 
elastic behavior. This function generally represents one or a com- 
bination of stress invariants and can be a smooth continuous func- 
tion (von Mises) or discontinuous piecewise linear (Tresca) . 

In the case of subsequent yielding from a plastic state, the 
function used to define the elastic limit is referred to as the 
subsequent yield or loading function, and the corresponding yield 
condition can be represented as 

f(c..,a..) - (a ) n = 0 , (B“l) 

ij lj o' v ' 
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where £ is a homogeneous function of order n of its arguments, 
o : j_j are the stress components, ctj_j represents a measure of the 
degree of work hardening, and a Q is the yield stress. 

The yield and loading conditions serve to establish criteria 
for loading, unloading, or neutral loading from elastic or plastic 
states, respectively (Ref. 25) . Additional information, in the 
form of a constitutive relation between increments of plastic 
strain, stress, and stress increment, is required to describe the 
plastic behavior of a material. This constitutive relation, termed 
the flow rule, is based on the maximum plastic work inequality 
(Ref. 69): 



a. . ) de . . >0 

ij - 


(B-2) 


where depj represents the plastic strain increment components 

resulting from the stress state opj on the subsequent loading 

surface,, and a?. is any other stress state on the surface be- 
ij 

neath it. In addition to providing the convexity condition on the 
loading surface, the inequality leads to the normality condition 
associated with the plastic strain increment vector. Thus, the 
flow rule is represented as 



dA 


5 f ( a ij' a ij> 

da. . 


(B-3) 


where dA is a positive scalar quantity. 

Having selected a yield condition and flow rule, we must now 
choose a function that will establish conditions for subsequent 
yielding from a plastic state. Choice of a hardening rule depends 
on the ease with which it can be applied in the chosen method of 
analysis as well as on its capability of representing the actual 
hardening behavior of structural materials. These requirements 
and the necessity of maintaining mathematical consistency with the 
yield function constitute the criteria for final choice of a harden 
ing rule. An appraisal of some of the available hardening rules 
f ollows . 


Kinematic Hardening 


The hardening behavior postulated in this theory assumes that 
during plastic deformation the loading surface translates as a 
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rigid body in stress space, maintaining the size, shape, and orien- 
tation of the yield surface. The primary aim of this theory, due 
to Prager (Refs. 26 and 27), is to provide a means of accounting 
for the Bauschinger effect. Such a capability is of particular im- 
portance in the strength prediction of aircraft structures, since 
severe loadings tend to be random, and failure may occur as the 
cumulative effect of a small number of such loadings . 


An illustration of kinematic hardening, as applied in conjunc- 
tion with the von Mises yield curve in the crq, 02 plane, is pro- 
vided in Fig. 45. The yield surface and loading surface are shown 
here for a shift of the stress state from point 1 to point 2. De- 
noting the translation of the center of the yield surface by a - , 
we may represent the loading function f in the form f(a^j _a ij)j 
the subsequent yield condition is given as 





(B-4) 


As a consequence of assuming a rigid translation of the load- 
ing surface, kinematic hardening predicts an ideal Bauschinger ef- 
fect for completely reversed loading conditions. That is, the 
magnitude of the increase of yield stress in one direction results 
in a decrease of yield stress of the same magnitude in the. reverse 
direction. Kadashevitch and Novozhilov (Ref. 70) have independently 
developed a hardening rule almost identical to Prager' s. In their 
theory, the total translation of the yield surface is regarded as 
associated with "internal microstresses" that remain in the body 
upon unloading. It is these internal microstresses that are con- 
sidered to be responsible for the Bauschinger effect. 

The kinematic hardening theory as set forth by Prager postu- 
lates that the increments of translation of the loading surface in 
nine dimensional stress space occur in the direction of the ex- 
terior normal to the surface at the instantaneous stress state. 

This geometric relation can be expressed analytically by 


da . . = cde . . , 
ij 


(B-5) 


where de-jj is the increment of plastic strain, which, according 
to the flow rule, Eq. (B”3), is in the direction of the exterior 
normal to the loading surface, and c is a parameter character- 
izing the hardening behavior of the material. 
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As indicated in Refs. 71 -73, inconsistencies arise when the 
theory is applied in various subspaces of stress, that is, when 
the symmetry of the stress tensor or the fact that some of the 
stress components are zero is taken into account in reducing the 
number of dimensions of the stress space. 

These inconsistencies produce the result that the loading sur- 
face will not, in general, translate in the direction of the ex- 
terior normal in a subspace of stress when it is made to do so in 
the full nine dimensional stress space. Reference 72 specifies 
stress conditions under which a linear transformation of variables 
enables the loading surface, in the transformed subspace, to trans- 
late in the direction of the exterior normal. 

To avoid this difficulty in implementing complete kinematic 
hardening, Ziegler (Ref. 28) has proposed a modification of Prager's 
rule. This replaces Eq. (B-5) with the following expression for 
the increment of translation 

da. . = dp (a. . “ a. .) dp > 0 . (B-6) 


The geometrical significance of this modification is shown 

(P) 

in Fig. 46, where the increment of translation da;. 7 , computed on 


the basis of Prager's rule, is compared with the increment of 

translation da , computed on the basis of Ziegler's modified 

rule. Note that in the latter case the increment of translation 

da£?) is in the direction of the vector from the center of the 
ij 

yield or loading surface to the point representing the instanta- 
neous stress state. 


The scalar dp in Eq. (B”6) is determined from the condition 
that the stress state must remain on the translated loading sur- 
face during plastic deformation. From Fig. 47 it can be seen that 
this condition may be represented as 



da. .) 

-i -i ' 


df 


ir da. . 
J 3-J 


= 0 


Substituting Eq. 


(B-6) into Eq. 


(B-7), we have 


(B-7) 
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d|_L = 


(B-8) 


df 

ag ii 


da. . 
ij 


(a. . - a. .) 
ij iJ 


df 

da. . 


An expression for the scalar factor dA associated with the 
flow rule, Eq. (B _ 3) , can be determined by recognizing that the 
vector cdc^j , shown in Fig. 47, is the projection of da^j on the 
exterior normal to the loading surface at the instantaneous stress 
state. Thus, we can write 


df 


(da.. - cds..) r = 0 
ij il da.. 


ij' 


ij 


(B-9) 


Using the flow rule, Eq. (B“3) , to substitute for de^j in 
Eq. (B“9) results in the following expression for dA 



(B-10) 


where the summation convention is adopted in 9-space. The flow 
rule now becomes 



(B-ii) 


Isotropic Hardening 

This theory assumes that during plastic flow the loading sur- 
face expands uniformly about the origin in stress space, maintain- 
ing the same shape, center location, and orientation as the yield 
surface. Figure 48 illustrates, on the basis of a simplification to 
a two dimensional plot, the yield and loading surfaces when the 
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stress state shifts from point 1 to 2 . Unloading and subsequent 
reloading in the reverse direction result in yielding at the stress 
state represented by point 3. The path 2-3 will be elastic, and 
0-2 is equal to 0-3. 

It can be seen that the isotropic type of representation of 
work hardening does not account for the Bauschinger effect exhibited 
by most structural materials. In fact, this theory provides that 
because of work hardening the material will exhibit an increase in 
the compressive yield stress equal to the increase in the tensile 
yield stress. Furthermore, since plastic deformation is an aniso- 
tropic process, a theory that predicts isotropy in the plastic 
range cannot be expected to lead to realistic results when non- 
proportional loading paths (not necessarily completely reversed) 
are considered. This conclusion has been indicated experimentally 
in Refs. 74 -77. For monotonic loading conditions, however, iso- 
tropic hardening is satisfactory and is commonly used. 

In the case of isotropic hardening, f(a-j_j) in Eq. (B“l) re- 
tains its original form, except that a Q increases in magnitude 
to provide for the expansion of the yield surface as work hardening 
proceeds. Equation (B - l) then serves to define a Q , which is re- 
ferred to as the "effective stress." The corresponding relation- 
ship between plastic strain increments and stress increments for 
the case of isotropic hardening is obtained by introducing the ex- 
pression for f(crpj) in Eq. (B - l) into the flow rule of Eq. (B - 3) , 
with c 0 treated as a constant. 

Work-Hardening Moduli 

A theory that combines kinematic and isotropic work hardening 
would require the yield condition of Eq. (B _ l) to be modified as 

follows. 


f(a ij ' a ij } " F(A) = 0 ' (B-12) 

where F(A) is a measure of the expansion of the yield surface in 
stress space, and, during plastic deformation, A is a monotoni- 
cally increasing scalar function. Equation (B _ 12) reduces to the 
case of kinematic hardening when F(A) = constant, and to iso- 
tropic hardening when ct^ = 0 and F(A) is monotonically in- 
creasing. 
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From experimental considerations concerned with the behavior 
of metals under cyclic loading, it has been observed that a steady 
cycle of alternating plastic flow is reached after a certain number 
of cycles (Ref. 78) . For some materials under a fixed strain (or 
stress) range, fully reversed cyclic tension-compression tests may 
tend to increase the stress (or strain) range toward a fixed steady- 
state value (work stiffening) . Conversely, in some metals initially 
cold worked into the plastic range or heat treated, the stress (or 
strain) range can be lowered to the same steady-state value during 
plastic cycling (work softening) . The hypotheses of isotropic and 
kinematic hardening cannot, in general, predict this observed phe- 
nomenon. In the uniaxial case shown in Fig. 49, for instance, for 
stress cycling between the limits ±a , the isotropic hardening 

D 

model predicts completely elastic behavior beyond the first half- 
cycle, whereas kinematic hardening predicts steady cyclic plastic 
straining beyond the first cycle — neither of which may represent 
the actual behavior. 

In proposing a work-hardening model to account for the be- 
havior of metals under cyclic loading conditions, Mroz (Refs. 79 
and 80) has introduced the notion of a field of work-hardening 
moduli and the variation of this field during the course of plas- 
tic deformation. In this proposed model, a stress-strain curve of 
an initially isotropic material is represented by n linear seg- 
ments of constant tangent plastic moduli, as shown in Fig. 50. in 
stress space, this approximation can be represented by n hyper- 
surfaces f D , fp, ..., f n , where f Q is the initial yield surface, 
and fp to f n define regions of constant work-hardening moduli. 
Figure 51 illustrates these hypersurfaces in the op “ 02 plane 
for an initially isotropic material. As seen from this figure, 
the surfaces f Q , f p , ... are similar and concentric, and for sim- 
plicity are schematically represented by a family of circles . If 
we consider proportional loading in the 02 direction, correspond- 
ing to a in Fig. 50, and if we assume that the surfaces can ex- 
perience a rigid translation without experiencing a change of size 
or orientation, then when the stress state reaches point A on 
Fig. 51a, the surface f Q will translate until it reaches the cir- 
cle fp at the stress corresponding to point B. The circles f Q 
and fp translate together until the point C is reached, where 
now ^o’ ^1’ anc * ^2 are attac hed at a common point of contact. 

For unloading and subsequent reversed loading, when the stress 
reaches a point corresponding to point E (Fig. 51b), reverse plas- 
tic flow occurs and the surface f Q translates downward along the 
o axis until it reaches the surface fp at F. Mroz further 
proposes that the curve of reverse loading in Fig. 50 joins the 
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curve OA / B ' G that is obtained by symmetry with respect to the 
origin from OABC. Thus, the curve of reverse loading EFG is 
uniquely defined by the curve of primary loading, represented by 
an equation of the form a = f(e). If a new coordinate system 
(c,e) with origin at C is used, we have for the curve CEFG 

I a = f (I e) . (B-13) 


In the generalization of this model to nonproportional load- 
ing, it is assumed that during translation of the hypersurfaces 
the individual surfaces do not intersect but consecutively contact 
and push each other. If we consider two neighboring surfaces fp. 
and having centers^a^ 0^ and Op + p defined by the posi- 

tion vectors a.. and a.. , then the equations of these two 
ij iJ 

surfaces are 



k k+1 

where a_ and a n are constants used to define the size of the 
o o k 

surface. From Ref. 80, the instantaneous translation of f ( cp- ^ ” a ij) 


is given by 
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where 


dp = 


df/dapjdapj 



On reaching surface fp^q, the position of fp- can be obtained 
from the following relation: 


k+1 

k+1 k+1 °o / k 

a. . - a. . = — r— a., 

ij il a k v 



(B-16) 
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i.e., the vectors connecting the center of the surface with the 
stress point are parallel for the surfaces fk and fk+1 • When 
the centers of the surfaces originally coincide, the mathematical 
description of the translation of the yield surface is identical 
to that given by Ziegler (Ref. 28) and presented here in Eq. (B“6) . 


The work-hardening modulus for each hypersurface is developed 
in Ref. 80 from a flow rule in the form of Eq. (B“3) , and is given 
as 


c 


jg ii d£ ii 

de . . de . . 
11 ij 


(B-17) 


for multiaxial stress. Equation (B - 17) represents a generalization 
of the plastic tangent modulus in the uniaxial stress state. A 
corresponding modulus of the form c = h(cr,e), where h is a 
functional representation of the instantaneous slope of an effec- 
tive stress (a) -effective strain (e) curve, is presented in 
Ref. 25. It should be noted that when f^ tends to infinity, so 
that the work-hardening modulus is constant (the work-hardening 
curve being represented by a straight line) , the theory proposed 
by Mroz is identical to Prager's kinematic hardening model. 

The further generalization of the theory of work-hardening 
moduli is associated with an expansion or contraction of the sur- 
faces f Q , f-^, ..., so that transitory phenomena (work stiffening, 
work softening or nonisothermal conditions) can be treated. Thus, 
the hypersurfaces fk are not constants but functions of a scalar 
parameter, A, monotonically increasing during plastic flow. One 
suggestion for A is presented in Ref. 80, where it is proposed 
that 


k 

c 

o 




y 


where 


A = 


0 


(de . .de . .) 2 dt , 
i] ij 


and t represents an interval of time or load, etc. 


(B-18) 
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The effects of cyclic loading, as predicted by kinematic 
hardening or the theory of Mroz, can be expected to give, at best, 
a simplified approximation of the actual behavior of structural 
metals under cyclic loading. These theories idealize the behavior 
of the location of the yield surface (s) in a general direction of 
prior plastic deformation. Experiments have shown that the be- 
havior of subsequent yield surfaces are far more complex than a 
mere translation and/or expansion of the original surface. How- 
ever, we believe that these theories satisfy our original criteria 
requiring mathematical simplicity and capability of predicting the 
essential features of cyclic plastic behavior. Results for sev- 
eral of those problems in which strain hardening is considered in 
this report were obtained by using the Prager-Ziegler kinematic 
hardening theory with a constant work-hardening modulus. For 
others, the concept of fields of work-hardening moduli was used 
where, for a point undergoing loading in the plastic range, a new 
modulus is computed for each increment of load. Further, for re- 
versed loading, a new coordinate system, used to define the re- 
versed stress-strain behavior, is chosen at a point corresponding 
to E in Fig. 51b. 
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COAL OF INVESTIGATION APPROACH &E5SZLX 



Fig. 1 Chronological Summary of "Nonlinear Finite Element 
Analysis Studies at Grumman 
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Type of Analysis 


Loading 


Application 


I. 

Membrane Stress 

® Arbitrary concentrated and/or 
distributed in-plane (membrane) 
loads . 

• Unloading from any plastic 
state and reversed loading. 

Thin walled membrane stressed 
structures, such as sheets with 
arbitrarily shaped holes and 
cutouts; general out-of -plane 
thin -walled membrane structures. 

II. 

1. Pending 

2. Combined bending and 
membrane 

• Arbitrary concentrated and/or 
distributed normal and in- 
plane loads. 

• Unloading from any plastic 
state and reversed loading. 

Plates of arbitrary shape with 
holes and cutouts; stiffened 
panel; general out-of -plane thin 
walled structures 

III. 

Axisymmetric bodies 
of revolution 

• Axisymmetric line or distributed 
loads . 

« Unloading from any plastic state 
and reversed loading. 

Thin walled axisymmetric stif- 
fened and unstiffened shells of 
revolution; thick walled axisym- 
metric structures; junctions be- 
tween branched structures; transi- 
tional shells (thin - thick). 

IV. 

Three dimensional 

« Arbitrary concentrated and/or 
undistributed loads . 

• Unloading from any plastic 
state . 

Behavior of arbitrarily shaped 
three dimensional bodies , 

V. 

Composite plates 
under plane stress 

e Arbitrary concentrated and/or 
distributed in-plane loads. 

• Unloading from any plastic 
state and reversed loading. 

Laminated composite plates under 
generalized plane stress condi- 
tions . 


Fig. 2a Nonlinear Analysis Capabilities of the PLANS System 

























-THEOR. (elastic) Ref (38) 

FE. ( elastic ) 

FE. (elasto-plastic )a 0 =.376 cr y j e j ^ 
EE. (elasto plastic ) cr Q = .753cr yie | d 
F.E. residual cr 0 = .376ayj e |d 
FE. residual a n = .753 a\. io M 


a o = •^^yield 5Ta. 

a 0 =.376oyj e |d 0 i 

a o ~ ‘’’yield 2 CM 



■feral 



1.0" 15“ (x-a) 0 


a/b = 0.096 
b/c = 1.300 



□ □ □ 


/a and 


(a) Stress Distribution Along y = 0 (b) Displacement Profile Along Crack Length 


•PLASTIC ANALYSIS OF A UNIFOPxMLY LOADED SHEET WITH A CENTRAL CF 









(a) PLASTIC ZONES flooding) (b) PLASTIC ZONES (unloading) 

Fig. 5 PROPAGATION OF PLASTIC ZONE IN A UNIFORMLY 
LOADED SHEET WITH A CENTRAL CRACK. 




Fig. 6 IDEALIZATION OF QUADRANT OF FULL SPHERE 
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w 0 /h 

Fig. 13 PRESSURE VERSUS NORMAL DISPLACEMENT, w , AT APEX ( ^ = 0°) 
FOR DIFFERENT LOAD INCREMENTS - TOR I S PHER I CAL SHELL 





00 

■C' 


10 



Fig. 16 LOAD VERSUS CENTER DEFLECTION OF A UNIFORMLY 
LOADED CLAMPED CIRCULAR PLATE 




(a ) Normal displacement at the hole boundary. 

Fig. 17 LOAD-DEFLECTION CURVES FOR RING -STIFFENED 
SPHERICAL SHELL UNDER EXTERNAL PRESSURE 
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(b) Normal displacement at r = 2.5 in =6.35cm 
Fig. 17 LOAD-DEFLECTION CURVES FOR RING-STIFFENED 
SPHERICAL SHELL UNDER EXTERNAL PRESSURE 
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- 4 . 



(a ) Distribution of circumferential stress resultant for unstiffened hole. 



0 .25 .50 .75 1.0 r/a 

(b) Distribution of circumferential stress resultant for stiffened hole. 


Fig. 18 RING-STIFFENED SPHERICAL SHELL UNDER EXTERNAL PRESSURE 
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( b ) Distribution of longitudinal stress 
Fig. 19 RECTANGUU 




y 


z 





ay along x axis , z = 0 (f = p/4bc ) 



.2 0 

; <r x along x axis, z= 0 (f=p/4bc } 


PRISM 




.6 


IN-ALUMINUM LAMINATE 


[ 90 / 90 / 0/0 ] 





Fig. 22 IDEALIZATION FOR PANEL WITH A CUTOUT 



Fig. 23 INTERLAMINAR SHEAR STRESSES AROUND CUTOUT 
IN BORON -EPOXY [+45] LAMINATES 



REF. 55 

O SOLUTION USING TANGENT MODULUS METHOD 
a SOLUTION USING EFFECTIVE LOAD METHOD 



Eh4 


FIG. 24. CENTRAL DEFLECTION VERSUS LOAD FOR A UNIFORMLY 
LOADED CLAMPED CIRCULAR PLATE 
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Eh 


SOLUTION USING TANGENT MODULUS METHOD 
SOLUTION USING EFFECTIVE LOAD METHOD 


BENDING 
STRESS AT 
EDGE 


LINEAR BENDING 

STRESS AT CENTER 
\ 



MEMBRANE 
STRESS AT 
EDGE 


w „ 


Fig. 25. STRESS VERSUS CENTRAL DEFLECTION FOR A 

UNIFORMLY LOADED CLAMPED CIRCULAR PLATE 


CENTRAL AXIAL LOAD, P (N) 



Fig. 26 LOAD VERSUS CENTRAL DEFLECTION FOR A 

CENTRALLY LOADED CLAMPED SPHERICAL CAP (X = 6) 
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LOAD 


ELASTIC, SMALL DEFL. LARGE DEFL & 



Fiq . 27 a LOAD VERSUS CENTER DEFLECTION FOR A SIMPLY- 
SUPPORTED CIRCULAR PLATE (2a/h =40.6) 






Fig. 27b CIRCUMFERENTIAL STRAIN DISTRIBUTION AT LOWER SURFACE 
FOR A SIMPLY-SUPPORTED CIRCULAR PLATE (2a/h =40.6) 





CIRCUMFERENTIAL MICROSTRAIN « 



RADIAL LOCATION R, cm 

Fig, 27c CIRCUMFERENTIAL STRAIN DISTRIBUTION AT UPPER SURFACE 
FOR A SIMPLY-SUPPORTED CIRCULAR PLATE (2a/h =40.6) 


97 






E = 30 X 10 psi = 20.6 X IO ,w N/m 1 " 
v = 0.3 

« = 0.1 
0 = .002 

r --,1/4 1/2 

X * 2 [3(1 - 1/ (H/h) 



FIG. 28. ELASTO- PLASTIC BUCKLING OF SIMPLY-! 
UNDER UNIFORM EXTERNAL PRESSURE (L< 
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CENTER DEFLECTION/ t 


29 LOAD VS. CENTER DEFLECTION OF A SIMPLY- 
SUPPORTED CIRCULAR PLATE 


Axial Stress Resultant, N x x 10 


R /t =120 
L = 25 in 



E la stic , Ideally plastic 

Elastic 

R/t =120 

E = 10 x I0 3 ksi = 6.9 x I0 10 N/m 2 
v = 0.3 

®"yield = 30 ksi =2.1 x I0 7 N/m 2 
L = 25 in =^0.64 m 



Co) Load vs axial end displacement. { b) E lastic and plastic rsdial displacement 

profile. 


Fig. 30 SIM PLY -SUPPORTED, STIFFEN ED CIRCULAR CYLINDER 
SUBJECTED TO COMPRESSIVE AXIAL LOADS. 
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(a) Applied end load vs end deflection for clamped and simply-supported circular cylinders 



(b) Radial shell displacement along axis for clamped and simply-supported end loaded cylinders. 
Fig .31 UNSTIFFENED CIRCULAR CYLINDERS 
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01 * ui / N 



fa) Axial load versus axial end defleclion. 

Fig. 32 TRUNCATED CONICAL SHELL 


« 17.090® 
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(b) Norma! she!! displacement along shell axis. 

Fig. 32 (cont.) CLAMPED TRUNCATED CONICAL SHELL 



a) Main Program b) ELAS 


Fig. 33 


Basic Flow of Each Analysis 


Program in "PLANS" 
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c) PLAS 


Fig. 33 (Cont . ) Basic Flow of Each Analysis Program in "PLANS" 






























Displacement Assumption : 
u(x,y)= a, + a 2 x+a 3 y 
v (x,y )= a 4 + a 5 x+a 6 y 

Initial Strain Distribution: 
(x,y) = constant 

x 


Fig. 35 CONSTANT STRESS TRIANGLE (CST) 



Fig. 36 LINEAR STRAIN TRIANGLE ( LST ) 




Fig. 37 WARPED SHEAR PANEL (GARVEY) 



(a) Constant strain element ( b) Linear slroin element 

Displacement Assumption^ 
u = a j + a 2 x u = a|+a 2 x+a 3 x 

Initial Strain Distribution 

e = constant 

Fig. 38 STRINGER ELEMENTS 
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X 



Displacement Assumption: 

w = a , + a 2 x t a 3 y + a 4 x 2 -t- a 5 xy + • • • + a 2) y 5 

2 3 

u = b| + bgx +b 3 y + b 4 x +b5xy+- +bioy 
v = C| +c 2 x tCjytc 4 x 2 + c 5 xy + "- + C| 0 y 3 
Initial Strain Distribution » 

€ (x,y,Z) = f(z) € j (z) 

i = I 

Fig. 39 TYPICAL TRIANGULAR ELASTIC - PLASTIC 
PLATE ELEMENT (PURE B) 


m xi 



X 

m zi 


f t f y' 

XI 1 


'hi 



m xj 


Displacement Assumption: 



2 




u 

V 

( 1) 

= 2^H oi u) 

Uj 

V i 

(1) 

t- H ti (e) * 

u.xi 

v >xi 

w 

i * 1 
2 

Wj 

w *xj 


W- 

i * I 


H = Hermitian polynomial of order k 
Fig. 40 BEAM ELEMENT 
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Displacement Assumption 


P = 


u 


1 

"0 

o 

o 

i 


_ 

u 

V 

= 

0 P 0 


V 

w 


0 0 P. 


w 


P l P 2 P 3- 


1 8 20 

Pj = — ( I + Cj C ) ( I * “*? j *? ) ( I *£j £ + *?j ’? + Cj C ~ 2 ) i = 1,8 
for a typical mid-edge node ( e-g. €; = 0, = 1 1, £ s = ± I ) 


P = V) ( 1+ t\t ) 


Plastic Strain Distribution 
€ j j = uniform 


Fig.43 ISOPARAMETRIC HEXAHEDRA 
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(a ) Laminate Model 


Layeri'- 


Layeri 



Fiber-Bearing 
xy i Orthotropic 
'/. Y.n| Segment 

, vz 

Interlaminar 
Shear 
Segment 


-M+i 


Layer/+I 

(b) Finite-Element Model 
Fig. 44 COMPOSITE ELEMENT 
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(p) 

da jj PRAGERS KINEMATIC HARDENING 
. (Z) 

da j j ZIEGLERS MODIFICATION 
a LOADING 



Fig. 45 KINEMATIC HARDENING Fig.46 COMPARISON OF PRAGER'S 

RULE WITH ZIEGLER'S MODIFICATION 


LOADING 




Fig. 47 HARDENING RULE AND FLOW Fig. 48 ISOTROPIC HARDENING 
LAW FOR WORK-HARDENING MATERIAL 
USING ZIEGLER'S MODIFICATION 
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